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I.  INTRODUCTION 

A.  DESCRIPTION  OF  THE  PROBLEM 

A  statistical  design  is  a  plan  according  to  which  an  experiment  is  patterned.  It 
provides  the  basis  upon  which  appropriate  statistical  tests  and  inferences  can  be  made 
after  the  experiment  has  been  performed.  The  selection  of  the  experimental  design  to 
be  used  in  a  given  situation  is  extremely  important  because  it  plays  a  predominant  role 
in  the  efficiency  of  the  experiment,  the  precision  with  which  the  objectives  are  met,  and 
the  total  effort  (and  cost)  expended  upon  the  experiment. 

The  concern  of  this  thesis  is  the  development  of  power  algorithms  for  several 
common  tests,  including  the  F-test  and  t-test.  The  experimenter  wants  a  test  to 
achieve  the  correct  decision  with  as  high  probability  as  possible  in  practical  operational 
tests  and  evaluations  given  certain  conditions  hold. 

B.  SCOPE  OF  THE  THESIS 

An  interactive  computer  program  is  developed  which  uses  power  algorithms  for 
several  common  tests.  The  program  can  be  used  by  decision  makers  who  may  not 
have  deep  knowledge  of  statistics  or  may  be  unfamiliar  with  the  details  of  the  design  of 
experiments.  The  output  helps  the  decision  maker  design  experiments  appropriate  for 
testing  hypotheses  with  common  statistical  tests.  The  experiment  designer  can  use  the 
program  to  evaluate  the  effect  upon  power  of  varying  underlying  design  parameters 
such  as  sample  size. 

Because  the  underlying  statistical  noncentral  distributions  used  in  computing 
power  do  not  have  tables,  approximation  methods  are  used  for  computing  the 
noncentral  CDFs.  Approximation  methods  were  developed  which  are  efficient, 
requiring  reduced  computer  CPU  time. 

C.  BACKGROUND 

A  hypothesis  is  a  statement  about  the  values  of  the  parameters  of  a  probability 
distribution.  For  example,  suppose  we  think  that  the  mean  yield  of  a  chemical  process 
is  more  than  94.5  percent.  Hypotheses  to  test  this  statement  might  be  expressed 
formally  as 

Ho  :  n  <  94.5 

Ha  :  ^  >  94.5 


The  statement  Ho  :  ^  :^  94.5  is  called  the  null  hypothesis,  and  Ha  :  [i  >  94.5  is  called 
the  alternative  hypothesis.  The  value  of  the  mean  specified  in  the  null  hypothesis 
might  typically  be  determined  in  one  of  two  ways.  It  may  be  the  result  of  some  theorv' 
or  model  regarding  the  process  under  study,  or  it  may  be  the  result  of  contractual 
specifications. 

To  test  a  hypothesis  one  usually  devises  a  procedure  for  taking  a  random  sample, 
computes  an  appropriate  test  statistic,  and  then  rejects  or  declines  to  reject  the  null 
hypothesis  Ho,  depending  on  the  outcome  on  the  test  statistic.  Part  of  this  procedure 
is  the  specification  of  the  set  of  values  for  the  test  statistic  which  lead  to  rejection  of 
Ho.  This  set  of  values  is  called  the  critical  region  or  rejection  region  for  the  test. 

Two  kinds  of  errors  may  be  committed  when  testing  hypotheses.    If  the  null 
hypothesis  is  rejected  when  it  is  true,  then  a  type  I  error  has  occurred.    If  the  null 
hypothesis  is  not  rejected  when  it  is  false,  then  a  type  II  error  has  been  made.    The 
probabilities  of  these  two  errors  are  given  special  symbols: 
a  =  P(type  I  error)  =  P{ reject  Ho  |  Ho  is  true) 
P  =  P(type  II  error)  =  P(fail  to  reject  Ho  |  Ho  is  false) 

Generally,  the  determination  of  P  requires  the  use  of  a  noncentral  sampling 
distribution  (e.g.,  noncentral  F,  noncentral  x  .  and  noncentral  t)  for  which  tables  are 
not  readily  available.  The  noncentral  sampling  distribution  depends  on  a  parameter 
called  the  noncentrality  parameter.  The  noncentrality  parameter  is  usually  a  measure 
of  the  distance  (in  some  sense)  between  the  values  of  the  parameter  under  the  null  and 
alternate  hypotheses.  Thus,  type  II  error  rates  depend  upon  parameters  of  the 
distribution  of  the  test  statistic  under  the  alternative  hypothesis.  The  power  of  the  test 
for  a  specified  value  of  the  hypothesized  parameter  is  the  probability  the  test  would 
reject  Ho  when  Ho  is  false,  or  1  —  p. 

In  addition,  whatever  the  test  procedure  is,  rejecting  Ho  when  Ho  is  not  true  is 
usually  something  we  want  a  test  to  achieve  with  as  high  a  probability  as  possible. 
Therefore  we  want  the  power  of  the  test,  for  a  given  value  of  non-centrality  parameter, 
to  be  high.  Power  algorithms  are  used  to  assist  in  the  "best"  selection  of  parameters 
for  a  statistical  design.  There  do  not  exist  simple  closed  form  equations  for  computing 
the  required  probabilities  from  the  common  noncentral  sampling  distributions. 
Consequently,  one  must  use  numerical  approximations  to  determine  power.  In  this 
thesis  w^e  use  various  approximation  methods  to  compute  the  cumulative  probabilities 
for  the  noncentral  F,  noncentral  x    ,  and  noncentral  t  distributions. 
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II.  NONCENTRAL  CHI-SQUARE,T,F  DISTRIBUTIONS 

A.       NONCENTRAL  CHI-SQUARE  DISTRIBUTION 

If  Wl,  W2,  ...,  Wn  are  independently  distributed  as  N{n-  ,1),  i=  1,2,  ...,n,  then 
T]  Wi  has  a  distribution  known  as  a  noncentral  /"  (n,  X  )  where  n  is  the  number  of 
degrees  of  freedom  and  ?^  =  X  M  i  is  the  noncentrality  parameter.  When  \^  \  ^  ^  2 
= ...  =  |H  j^  =  0,  then  X  =  0,  and  the  noncentral  x  (n,0)  reduces  to  the  usual  central 
X^  (n)  with  n  degrees  of  freedom.   The  cumulative  distribution  function  of  x"  (n,X  )  is, 


'p■^'^;1,^)  =R-[X-''c«,x)  $  x  ]  (2.1) 


while   F(x;n,  X)  =  0   for  x  <  0 

It  is  possible  to  express  F{x;n,X)  for  x>0,  in  an  easily  remembered  form  as  a 
weighted  sum  of  central  x  CDF's  with  weights  equal  to  the  probabilities  of  a  Poisson 
distribution  with  expected  value  X/2.   That  is, 


ao 


5=0 


♦  p-l  A-CM+ij)   ^  ^J  (2.2) 


7  7 

Thus  a  X    (n,  X)  variable  can  be  regarded  as  a  mixture  of  central  x    variables.    This 

interpretation  is  often  useful  in  deriving  the  distribution  of  functions  of  noncentral  x 

random  variables. 

The  probability  density  function  can,  similarly,  be  expressed  as  a  mixture  of 

central  x    probability  density  functions 


A^x^'^  ,        .      e''^x- 


.2. 


JL 


j^o  r^^^^^^rf^j)  z't 


(2.3) 
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The  mean  and  variance  of  the  distribution  are 

E(X)=n+2>.   and   Var(X)=  2n+ 8^  [Ref.  l:p.  130]. 

B.       NONCENTRAL  F  DISTRIBUTION 

If  Yj  and  Y2  are  independent  and  Yj  is  x  j     (Hj  ,'k  )  and  Y2  is  X     (rij  )  then 

V=(Yi/n^)/(Y2/n2)  (2.4) 

is  distributed  as  F'(nj  ,n2  ,X  ),  the  noncentral  F  distribution  with  n^  and  n^   degrees  of 
freedom  and  noncentrality  parameter  X  .   Its  density  function  is  given  by 

•f(v;=  CLl ■ :  v^^o        ,. ,, 

J  K=o  p  (^ii-Kj  ir\^  +  ^\r)*<  >  ^2.5) 

where  -  >    -2]      4^      >it 


p(-^)(ivtyi,ir;^-i-^ 


when   X  =0  this  reduces  to  the  density  function  of  the  central  F(n,,n2  )  distribution. 
The  mean  and  variance  of  the  distribution  are 


E(v)=~^l  +  -±^'\         forn  >  2 


II  t- 

and 


Var(V)  =  ^ 


MVl^a^^jVM^ 


Lcnz-i)(Ti^-+)     »^- "^ 


■) 


(for  n^  >  4) 


When     X    =0,    these    reduce    to    the    mean    and    variance    of  the    central    F(npn2) 
distribution.   Derivation  of  (2.5)  is  shown  in  [Ref  l:p.  189]. 

The  cumulative  distribution  of  V  can  be  expressed  in  terms  of  an  infinite  series  of 
multiples  of  incomplete  beta  functions,  as  follows: 


P^  ^  ^  ^^  ^'^  %S  ( ^  e-')  •  I^,  (i"'%  ^) 
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where   I  yi,       (— rtitj    —  )        is  the  incomplete  beta  function. 

In  this  thesis,  we  use  Paulson's  Approximation  from  Severo  and  Zelen  [Ref  2] 
Their  approximation  is  as  follows: 

Pr(V  <  fp)  ~  a)(x)  (2.7) 


1 


where  >,  x      I    _  ^ 

X  =. ,    ,  and 

<I>  is  the  standard  Normal  CDF. 

The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  Pr(V  <  fo  ).  It  can  be  seen  in  TABLE  1  that  this 
approximation  gives  about  2  decimal  place  accuracy.  Thus,  this  approximation 
provides  adequate  accuracy  for  computing  the  power  of  F-tests  [Ref  3:p.  202]. 

The  program  we  used  in  the  computation  of  noncentral  F  CDF  values  is  shown 
in  Appendix  D. 

C.       NONCENTRAL  T  DISTRIBUTION 

The  ratio 

T  =  (U  +  ).)/V(Y/n)  (2.8) 

where  U  and  Y  are  independent  random  variables  distributed  as  standard  normal 
(N(0,1))  and  x     with  n  degrees  of  freedom  , respectively,  is  said  to  have  a  noncentral  t 
distribution  with  n  degrees  of  freedom  and  noncentrality  parameter  X.     Sometimes   X, 
(or  even  .5X  ),  rather  than  X,  is  termed  the  noncentrality  parameter.    If  X  is  equal  to 
zero,  the  distribution  is  central  t  with  n  degrees  of  freedom. 
The  cumulative  distribution  of  T  is  [Ref  l:p.  201] 
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» 


Pr(T  <  I)  =- 


^t- 1  n  a 


r© 


X      e       F=r=^        e     ^       ^^-'Wx       (2.9) 


TABLE  1 
CDF  OF  A  NONCENTRAL  F  DISTRIBUTION 


nl 

n2 

X 

fo 

approx 

exact 

3 

10 

4 
4 
16 

16 

3.708 
6.552 

3.708 
6.552 

0.7499 
0.9192 
0.2023 
0.5190 

0.745 
0.918 
0.206 
0.517 

3 

20 

4 
4 
16 
16 

3.098 
4.938 
3.098 
4.938 

0.7067 
0.8894 
0.1186 
0.3488 

0.700 
0.887 
0.126 

0.347 

5 

10 

6 
6 

24 
24 

3.326 
5.636 
3.326 
5.636 

0.7337 
0.9143 
0.1553 
0.4629 

0.731 
0.914 
0.158 
0.461 

5 

20 

6 
6 

24 
24 

2.711 

4.103 
2.711 
4.103 

0.6685 
0.8715 
0.0643 
0.2437 

0.664 
0.870 
0.069 

0.245 

8 

10 

9 
9 

36 
36 

3.072 
5.057 
3.072 
5.057 

0.7159 
0.9087 
0.1166 

0.4088 

0.714 
0.908 
0.119 
0.408 

8 

30 

9 
9 

36 
36 

2.266 
3.173 
2.266 

3.173 

0.581 

0.8157 
0.0146 
0.0846 

0.578 
0.813 
0.017 
0.088 

In  this  thesis,  we  use  an  approximation  method  suggested  by  Abamowitz  and 
Stegun  [Ref.  4:p.  949]. 


Pr(T  <  t)  -  CD  (x) 


(2.10) 


where  x  = 


t(i-^)-x 


[^^^)" 


and  <I)  is  the  standard  Normal  CDF. 
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The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  Pr{T  <  to  ).  It  can  be  seen  in  Table  2  that  this 
approximation  gives  about  2  decimal  place  accuracy.  Thus,  this  approximation 
provides  adequate  accuracy  for  computing  the  power  of  t-tests. 

The  program  we  used  in  the  computation  of  noncentral  t  CDF  values  is  shown 
Appendix  E. 

TABLE  2 
CDF  OF  A  NONCENTRAL  T  DISTRIBUTION 


n 

X 

to 

approx 

exact 

4 

0.112 

7.944 

2.24 

18.494 

4 
4 
10 
10 

0.017 
0.993 
0.026 
0.993 

0.01 
0.99 
0.01 
0.99 

24 

5.786 

13.981 

16.118 

33.078 

9.797 
9.797 
24.494 
24.494 

0.012 
0.991 
0.014 
0.992 

0.01 
0.99 
O.Ol 
0.99 
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III.  POWER  OF  THE  F-TEST 

A.  INTRODUCTION 

Under  appropriate  conditions,  the  best  test  for  testing  equality  of  several  means 
is  the  analysis  of  variance  test.  Analysis  of  variance  has  a  wide  application.  It  is  one 
of  the  most  useful  techniques  in  the  field  of  statistical  inference.  As  in  any 
hypothesis-testing  situation,  the  power  of  the  F  test  is  of  interest  to  the  experimenter. 

In  this  chapter  we  will  discuss  the  power  of  F  tests  and  provide  an  example.  To 
give  an  overall  evaluation  of  the  power  of  F  tests  in  the  analysis  of  variance,  we  may 
use  power  curves.  An  important  use  of  the  power  curve  is  to  guide  the  experimenter  in 
selecting  ^he  sample  size  (number  of  replicates)  so  that  the  design  will  be  sufficiently 
sensitive  to  important  potential  differences  in  the  treatments.  We  will  consider  power 
curves  for  one-way  and  multi-way  analyses  of  variance  (ANOVA's). 

B.  THE  ONE-WAY  CLASSIFICATION  ANALYSIS  OF  VARIANCE 

Suppose  we  have  k  different  levels  of  a  single  factor  that  we  wish  to  compare. 
The  different  levels  of  the  factor  are  often  called  treatments.  The  observed  responses 
from  each  of  k  treatments  is  a  random  sample  on  a  random  variable.  The  data  would 
appear  as  in  Figure  3.1. 


Treatment 

1 
2 

k 

Observation 

Yll 

Y12     ... 

hn 

Y2I 

Y22 

\n 

. 

hi 

hz      ■■■ 

\n 

Figure  3.1     Typical  Data  for  One-way  Classification  ANOVA. 
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We  will  find  it  useful  to  describe  the  observations  by  the  linear  statistical  model: 

Yj-  =  H  +  T-  +  cj-  (3.1) 

where  Y-  is  the  (ij)      observation,  ^  is  a  parameter  common  to  all  treatments  called 
the  overall  mean,  T-  is  a  parameter  peculiar  to  the  i      treatment  called  the  i      treatment 
effect,  and  z-  is  a  random  error  component,  assumed  to  be  IID  N  (0,  <t    ). 
The  indices  used  are: 

•  i  =  the  treatments,  i  =   1,2,. ..,k 

•  j  =    the  replication  per  treatment,  j  =   l,2,...,n 

The  objective  in  the  ANOVA  is  to  test  appropriate  hypotheses  about  the 
treatment  effects.  The  variance  <7~  is  assumed  constant  for  all  levels  of  the  factor. 
This  model  is  called  the  one-way  classification  analysis  of  variance  because  only  one 
factor  is  investigated. 

We  are  interested  in  testing  the  equality  of  the  k  treatment  effects,  so  the 
appropriate  hypotheses  are 

Ho  :  Tj  =  ^2  =  ...  ^  "^k  ^  '^ 
Ha  :      '^i  ^  "^j   ^oi"  some  i,j  . 

That  is,  if  the  null  hypothesis  is  true,  then  each  observation  is  made  up  of  the  mean  yi 
+  T  plus  a  realization  of  the  random  error  £-. 

The  results  of  the  ANOVA  procedure  is  summarized  in  Table  3. 

We  may  use  the  expected  values  of  the  mean  squares  to  verify  that  Fo  (in  Table 
3)  is  an  appropriate  test  statistic  for  Ho.  From  the  expected  mean  square  we  see  that, 
in  general,  MSE  is  an  unbiased  estimator  of  <y  .  Also,  under  the  null  hypothesis,  MSt 
is  an  unbiased  estimator  of  a  .    However,  if  the  null  hypothesis  is  false,  then  the 

expected  value  of  MSt  is  (T^  4-(n  ^  (T-)^)/(k-l).    The  expected  value  of  the  mean 

2 
square  error  is  c^. 

Therefore,  under  the  alternate  hypothesis  the  expected  value  of  the  numerator  of 

the  test  statistic  (Fo)  is  greater  than  the  expected  value  of  the  denominator  and  we 

would  reject  Ho  with  values  of  the  test  statistic  which  are  too  large.    That  is,  we  would 

reject  Ho  if 

^o  "^  Pa'(k-l)'(N-k) 
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TABLE  3 
THE  ANALYSIS  OF  VARIANCE  FOR  A  ONE-WAY  CLASSIFICATION 


Source  of 
Variation 

Sum  of 
Squares 

Degrees 
Freedom 

of 

Mean 
Square 

Fo 

Between 
Treatments 

SSt 

k-1 

MSt 

MSt/MSE 

Error(within 
treatments) 

SSE 

k(n-l) 

MSE 

Total 

SST 

kn-1 

where  the  numerator  and  denominator  df  are  k-1  and  k(n-l)  =  N-k,  and  a  is  the  type  I 
error  rate. 

The  power  of  the  test  is: 

1-P  =  P(Fo  >  F„,(i^.j),(^.i^)  I  Ho  is  false)  (3.2) 

To  evaluate  the  P  in  Equation  3.2  we  need  to  know  the  distribution  of  the  test  statistic 
Fq  if  the  null  hypothesis  is  false.  It  can  be  shown  that,  if  Ho  is  false,  the  statistic  F^ 
has  the  noncentral  F  distribution  with  k-1  and  k(n-l)  degrees  of  freedom  and 
noncentraUty  parameter  X,  given  by 

I  =       

The  noncentraUty  can  be  interpreted  as  the  squared  standardized  distance 
between  the  origin  and  (Tpt2,  ...  ,Tj^).  The  ratio  ^  (t-  )  /  <y  is  called  the  squared 
standardized  distance.  If  only  an  estimate  of  a  is  available,  one  may  replace  <J  with 
the  estimate  [Ref.  5:p.  34]. 

C.       THE  MULTI-WAY  CLASSIFICATION  ANALYSIS  OF  VARIANCE 

Many  experiments  require  a  study  of  the  effects  of  two  or  more  factors.  It  can 
be  shown  that,  under  certain  conditions,  factorial  arrangements  are  the  most  efficient 
designs  for  this  type  of  analysis. 


One  of  the  simplest  factorial  experiments  involves  only  two  factors  or  sets  of 
treatments  say  factor  A  and  factor  B.  Suppose  there  are  a  levels  of  factors  A  and  b 
levels  of  factor  B,  and  these  are  arranged  in  a  2-way  factorial  design;  that  is,  each 
replication  of  the  experiment  contains  all  ab  treatment  combinations.  Assume  there 
are  n  replications  of  the  experiment,  and  let  Y;-u  represent  the  observation  taken  under 
the  i      level  of  factor  A  and  the  j      level  of  factor  B  in  the  k      replication. 

The  data  can  be  summarized  as  shown  in  Figure  3.2.  The  order  in  which  the  abn 
observations  are  taken  is  selected  at  random. 


Factor  B 


Factor 
A 


1      1       1        2     1   ... 

1   '^lll'^112   1   '^121''^122   1 
1  1  ,..•  ,  Yi2n    1   '•••  '  ^12n   I 

1   Y212.Y212    1   Y22i,Y222   1 
2  1   ,...  ,  Y2in   1   >•••  .   Y22n  1 

.  1    ...      1    ...      1   . 

1   ^all'^al2    1   ^a21'^a22   1 
a  1   ,.  . .  ,  Yg2n   '   '•  •  •  >  '''a2n   ' 

b 

^Ibl 

,  .  .  .  , 

'Ylb2  1 
"^Ibn    1 

'^2bl 
, . . .  , 

''^2b2  1 
^2bn    1 

•  • 

i 

^abl' 
, . . .  , 

^ab2  1 
^abn    1 

Figure  3.2    Typical  Data  notation  for  a  Two-way  Classification. 
The  observations  may  be  described  by  the  linear  model: 


Yijk=  ^^  +  ^i  +  Pj4-'^Pij+^ijk' 


(3.3) 


where 
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• 


\i  =  the  overall  mean  efTect; 


•  T  J  =  the  true  effect  of  i       level  of  factor  A,  i  =   1,2,. ..,a; 

•  P  :  =  the  true  efTect  of  j  ^^  level  of  factor  B,  j  =  l,2,...,b; 

•  T  P  •;  =  the  efTect  of  the  interaction  between  x  •  and  P  ;  ;  and 

•  c  iji^  =  a  random  error  component,  assumed  to  be  IID  N  (0,cr     ). 

Both  factors  are  assumed  to  be  fixed.  It  is  usually  assumed  that  the  treatment  effects 
are  defined  as  deviations  from  the  overall  mean,  so  ^  t  :  =  0  and  5^  P  ;  =  0- 
Similarly,  the  interaction  efTects  are  fixed  and  usally  defined  so  that  T]  (  '^  P  )  j;  =  0, 
where  the  summation  is  over  either  i  or  j.  Since  there  are  n  replicates  of  the 
experiment,  there  are  a  total  of  abn  observations. 

We  are  interested  in  testing  various  hypothesis  about  the  parameters  in  equation 
3.3  .  An  appropriate  hypothesis  testing  procedure  would  again  be  analysis  of  variance. 
More  specifically,  as  we  are  considering  two  controllable  sources  of  variation  (A  and 
B),  the  procedure  is  called  the  two-way  classification  analysis  variance. 

In  order  to  test  the  hypothesis  Ho  ;  T  •  =  0  for  i=  1,2,. ..,a  (no  row  factor 
effects).  Ho  ;  P  :  =  0  for  j=  l,2,...,b  (no  column  factor  efTects),  and  Ho  ;  (  T  P  )  •;  =  0 
(no  interaction  efTects),  we  can  express  the  total  sum  of  square  as: 

SST  =    SSA  +  SSB  +  SSAB  +  SSE.  (3.4) 

Here,  SSA  is  a  sum  of  squares  due  to  "rows"  or  factor  A,  SSB  is  a  sum  of  squares  due 
to  "columns"  or  factor  B,  SSAB  is  a  sum  of  squares  due  to  the  interaction  between  A 
and  B,  and  SSE  is  a  sum  of  squares  due  to  error.  The  degrees  of  freedom  associated 
with  each  sum  of  squares  are  shown  in  Figure  3.3. 

If  we  assume  e  -j^  are  IID  N(0,(T'^  )  and  apply  Cochran's  theorem  (Theorem  3-1) 

under  the  null  hypothesis  of  no  effects,  each  sum  of  squares  on  the  right-hand  side  of 

7  7 

Table  4  when  divided  by  d"^  is  distributed  as  X      with  the  indicated  number  of  degrees 

of  freedom,  and  these  statistics  are  independent. 

Theorem  3-1  (COCHRAN).  Let  Z  •  be  IID  N(0,1)  for  i=  1,2,... ,v  and  suppose 
^Z  i  Qj  +  Q2  "^  Q3  "•"•••  ■'"Qs  where  s  :^  v,  and  the  quadratic  form  Q-  has 

V-  degrees  of  freedom  (i=  l,2,...,s).  Then  Qj  Q2  ,..-,  Q3  are  independent 
chi-square  random  variable  with  Vi  V2  ..,  v  degrees  of  freedom,  respectively,  if 
and  only  if 

V  =  vj  +  V2  +  ...+  Vg. 
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Effect 

Degrees  of  freedom 

A 
B 

AB  interaction 
Error 

a-i 
b-1 
(a-l)(b-l) 
ab(n-l) 

Total 

abn-1 

Figure  3.3     Table  of  degrees  of  freedom  with  sums  of  squares.. 

Assuming  that  factors  A  and  B  are  fixed,  the  expected  values  of  the 
mean  squares  are: 


E{MSA)  =  (7^  +(  bn  5]  T  i  ^  )/(a-l)  ; 


(3.5) 


E(MSB)  =  (t2  +(bnXPi^)/{b-l); 


(3.6) 


E(MSAB)  =  (7^  +  (n  5:  X  (^  P  )  ii  ^  )/  (a-l)(b-l)  ;   and 


(3.7) 


E(iMSE)  =  (7^  . 


(3.8) 


Therefore,  to  test  the  hypothesis  Ho  :  T  j  =  0  ;  i=  l,2,...,a  (no  row  factor  effects),  and 
Ho  ;  p  ;  =  0  ;  j=  l,2,...,b  (  no  column  factor  effects),  and  Ho  :  (  i  P  )  ■:  =  0  (no 
interaction  effects),  we  would  divide  the  corresponding  mean  square  by  the  mean 
square  error.  Under  the  null  hypothesis  of  no  effect,  this  ratio  will  follow  an  F 
distribution  with  appropriate  numerator  degrees  of  freedom  and  ab(n-I)  denominator 
degrees  of  freedom,  and  the  critical  region  will  be  located  in  the  upper  tail.  The  test 
procedure  is  usually  summarized  in  an  analysis  of  variance  table,  such  as  shown  in 
Table  4. 
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TABLE  4 
THE  ANALYSIS  OF  VARIANCE  FOR  A  TWO-WAY  CLASSIFICATION 


Source  of 
Variation 

Sum  of 
Squares 

Degrees  of 
Freedom 

Mean 
Square 

Fo 

A  treatments 

SSA 

a-1 

MSA 

MSA/MSE 

B  treatments 

SSB 

b-1 

MSB 

MSB/MSE 

Interaction 

SSAB 

(a-l)(b-l) 

MSAB 

MSAB/MSE 

Error 

SSE 

ab(n-l) 

MSE 

Total 

SST 

abn-1 

To  the  compute  the  power  of  tests  in  two-way  ANOVAs.  the  procedure  is  the 
same  as  in  the  one-way  case.  Relationships  among  X,  the  numerator  degrees  of 
freedom  and  the  denominator  degrees  of  freedom  are  shown  in  Table  5  . 

In  a  similar  way  one  can  expand  to  the  general  multi-way  ANOVA  procedure 
[Ref  5:p.  124]. 

D.       THE  ALGORITHMS  AND  FLOWCHART 

A  program  for  computing  a  power  table  and  power  curve  in  general  ANOVA's 
(fixed  model)  is  shown  in  Appendix  F.  The  power  of  one-way  ANOVA's  are  solved 
with  these  programs,  using  the  following  the  sequence: 

1)  Determine  what  variable  to  include,  such  as  sample  size  vs  power,  number  of 
treatments  vs  power,  a-level  vs  power,  or  noncentrality  vs  power. 

2)  Detennine  the  values  of  the  maximum,  minimum,  and  increment  of  the  variable 
chosen. 

3)  Given  input  data,  compute  the  critical  value,  using  the  inverse  of  the  central 
F-distribution  (Appendix  A). 

4)  Compute  the  power  value  in  each  case;  that  is,  the  CDF  of  the  noncentral 
F-distribution  (Chapter  II). 

5)  Print  the  power  table  or  the  power  curve. 
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TABLE  5 
NONCENTRALITY  PARAMETERS  FOR  POWER  IN  A  TWO-WAY  ANOVA 


Factor  X 

Numerator 
DOF 

Denominator 
DOF 

A     — ^-T — 

a-1 

ab(n-l) 

b-1 

ab(n-l) 

AB  "l^^^j' 

(a-l)(b-l) 

ab(n-l) 

Multi-way  ANOVAs  have  the  same  algorithms  but  multi-way  ANOVA  may 
include  tests  of  interaction  efiects.  The  program  considers  only  up  to  three-way 
interaction  effects.  We  can  explain  how  to  compute  degrees  of  freedom  of  error  term 
in  m-way  ANOVA's,  assuming  the  balanced  case,  as  follows:  Total  degrees  of  freedom 
is  DOF(Total)=n  2^  k.(i)  -1,  where  k(i)  is  the  number  of  levels  of  the  i     factor  [Ref  9]. 

(1)  If  the  model  has  only  main  effects  without  any  interaction  effects, 

DOFl(Error)=  DOF{Total)  -  J]  (k(i)-l). 

(2)  If  the  model  has  several  factors  and  only  2- way  interaction  effects, 

DOF2(Error)=  DOFl(Error)  -  Y.  ^  (k(i)-l)(k(i)-l). 

(3)  IF  the  model  has  several  factors  and  only  3-way  interaction  effects, 

DOF3(Error)=  DOFl(Error)  -  SH  (k(i)-l)(k(j)-l)(k(k)-l). 

(4)  If  the  model  has  several  factors  and  2-way  and  3-way  interaction 

effects, 

DOF4(Error)  =  DOF3(Error)  -  ^  S  (k(i)-l)(k(i)-l). 
If  the  model  has  more  than  3-way  interaction  effects,  then  the  user  must  modify 
the  degrees  of  freedom  of  the  error  terms  accordingly.   A  flowchart  is  shown  in  Figure 
3.4. 
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Q         3TAPT ^ 


-> 


Deteririne  one-way  Of 
multi-way  ANOV  ? 


one-way 


Deternnne  what   to    include? 
.    sample   size   vs    power 
.    alpha-level    vs    power 
.    noncentral i ty  vs    power 


Compute  the  critical  value 
and  power 


Plot  power  curve 


_± 


Analyze  again? 


("         STOP         J 


m;uI  f  i  -way 


How  many  ways    (M)? 

For   i-factor   level  , 

How  many    level     (MFW)  ? 


Determine   how  to 
many   wo y    i rt eract inn?y 


Compute    CdrrtSLondi ny   to 
deqrees   of    treedoii 


Deterir.ine  what  to    include? 
.    •  of   replicates    vs    power 
.    noncent ral 1 ty  vs    power 
.    aloha-levtl    vs    power 


I 

/ 

Compute  the  critical 

value  and 

power 

Figure  3.4     System  flowchart  for  F-test  power. 
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E.       EXAMPLES  OF  F-TESTS 

1 .  Example  of  One-way  ANOVA  test  (number  of  replicates  vs  power) 

a.  Scenario 

A  manufacturer  suspects  that  the  batches  of  raw  material  furnished  by  his 
suppher  difier  significantly  in  calcium  content.  There  are  a  large  number  of  batches 
currently  in  the  warehouse.  Five  of  those  are  randomly  selected  for  study.  A  chemist 
wants  to  know  the  appropriate  sample  size  per  each  batch,  in  order  to  test 

Ho:  Tj  =t2  =  ...  =Ti^ 
against 

Ha  :      Tj  #  T;   for  some  i,j 

The  chemist  would  like  to  know  how  many  replicates  to  run  if  it  is  important  to  reject 
Ho  with  probability  at  least  0.9  when  the  standardized  distance  square  (  ^  t  •  /  <y  ) 
is  2  and  a  =  .05.  Thus  he  would  like  to  know  what  the  power  of  the  F-test  is  for  a 
range  of  possible  replicates.  He  decides  to  check  replicates  from  2  to  12  in  increments 
of  1. 

b.  Inputs 

1)  Select  one-way  ANOVA  test  (the  number  of  replicates  vs  power). 

2)  The  number  of  treatment  =  5. 

3)  a  -level  =  .05 

4)  Standardized  distance  square  =  2. 

5)  M.-^ximum  replicates  is  12,  minimum  replicates  is  2,  increment  is  1. 

c.  Start  program 

Do  you  want  to  analyze  one  way  ANOVA(y/n)? 

y 

Do  you  want  to  plot  n  (#  of  observation  per  treatment  vs  power)(y/n)? 

y 

Number  of  k(  #  of  treatment)? 

5 

Alpha-  level? 

? 

.05 
Standardized  squared  distance  value  ? 

? 
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Maximum  n  value  (the  maximum  value  on  X-axis)? 

7 

12 

Minimum  n  value  (the  minimum  value  on  X-axis)? 
(condition  :N  must  be  more  than  2) 

2 

Increment  n  value  ? 

1 

d.   Output 

The  output  is  presented  in  Table  6 

TABLE  6 
OUTPUT  OF  THE  ONE-WAY  ANOVA  EXAMPLE  (REPLICATES  VS  POWER) 


DOFl 

DOF2 

a 

F-INVERSE 

X 

POWER 

4 

5 

0.05 

5.29087 

4.00 

0.15373 

4 

10 

0.05 

3.52496 

6.00 

0.30024 

4 

15 

0.05 

3.10397 

8.00 

0.44935 

4 

20 

0.05 

2.91676 

10.00 

0.58677 

4 

25 

0.05 

2.81108 

12.00 

0.70332 

4 

30 

0.05 

2.74323 

14.00 

0.79562 

4 

35 

0.05 

2.69600 

16.00 

0.86449 

4 

40 

0.05 

2.66122 

18.00 

0.91328 

4 

45 

0.05 

2.63455 

20.00 

0.94631 

4 

50 

0.05 

2.61345 

22.00 

0.96776 

4 

55 

0.05 

2.59634 

24.00 

0.98118 

#  of  replicate 

POWER 

NN=     2. 

POWER  = 

0.15373 

NN=     3. 

POWER  = 

0.30024 

NN=     4. 

POWER  = 

0.44935 

NN=     5. 

POWER  = 

0.58677 

NN=     6. 

POWER  = 

0.70332 
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0.98118 
POWER 


0.56746 


0.15373 


NN  = 

7. 

POWER  = 

0.79562 

NN  = 

8. 

POWER  = 

0.86449 

NN  = 

9. 

POWER  = 

0.91328 

NN  = 

10. 

POWER  = 

0.94631 

NN  = 

11. 

POWER  = 

0.96776 

NN  = 

12. 

POWER  = 

0.98118 

2.000 

X-SCALE 


■+- 
7.000 
=  0.125E  +  00  UNITS 


+ + +-+ 

#  OF  REPLICATES     12.00 


Y-SCALE:   "\"=  0.138E-01  UNITS 
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2.  Example  of  one-way  ANOVA  test  (noncentrality  vs  power) 

a.  Scenario 

Five  brands  of  batteries  are  under  study.  It  is  suspected  that  the  life  (in 
weeks)  of  the  five  brands  is  different.  Five  batteries  of  each  brand  are  tested.  A 
manufacturer  wants  to  know  the  power  as  a  function  of  the  standardized  squared 
distances.  That  is,  would  like  to  know  what  the  power  of  the  F-test  is  for  a  range  of 
possible  standardized  distances  square.  He  decides  to  check  standardized  squared 
distances  from  I  to  5  in   .2   increment,  where  n=5  and  a  =  .05. 

b.  Inputs 

1)  Select  one-way  ANOVA  test  (the  standardized  distance  square  vs  power). 

2)  Number  of  treatment  =  5. 

3)  Number  of  replicates  =  5. 

4)  a  -level  =  .05. 

5)  Maximum  standardized  distance  square  is  5  and  minimum  standardized  distance 
square  is  1  and  the  increment  is  .2. 

c.  Start  program 
Start  program: 

Do  you  want  to  analyze  one-way  ANOVA(y/n)? 

y 

Do  you  want  to  plot  n(#  of  observation  per  treatment)  vs  power(y/n)? 

n 

Do  you  want  to  plot  alpha-level  vs  power(y/n)? 

n 

Do  you  want  to  plot  noncentrality  vs  power(y/n)? 

y 

#  of  observations  per  treatment  (n=  )  ? 

5 

#  of  treatment  (k=  )? 

9 

5 
Alpha-level  ? 

.05 
Maximum  standardized  squared  distance  value  range  ? 
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7 

5 

Minimum  standardized  squared  distance  value  range  ? 

7 
1. 

Increment  standardized  squared  distance  value  range  ? 

7 

.2 

d.   Output 

The  screen  output  (  Table  7  )  is  as  follows: 

TABLE  7 

OUTPUT  OF  THE  ONE-WAY  ANOVA  EXAMPLE  (NONCENTRALITY  VS 

POWER) 


DOFl 

D0F2 

a 

F-INVERSE 

\ 

POWER 

4 

20 

0.05 

2.91676 

5.00 

0.30382 

4 

20 

0.05 

2.91676 

6.00 

0.36314 

4 

20 

0.05 

2.91676 

7.00 

0.42203 

4 

20 

0.05 

2.91676 

8.00 

0.47946 

4 

20 

0.05 

2.91676 

9.00 

0.53461 

4 

20 

0.05 

2.91676 

10.00 

0.58677 

4 

20 

0.05 

2.91676 

11.00 

0.63550 

4 

20 

0.05 

2.91676 

12.00 

0.68050 

4 

20 

0.05 

2.91676 

13.00 

0.72163 

4 

20 

0.05 

2.91676 

14.00 

0.75885 

4 

20 

0.05 

2.91676 

15.00 

0.79223 

4  • 

20 

0.05 

2.91676 

16.00 

0.82192 

4 

20 

0.05 

2.91676 

17.00 

0.84812 

4 

20 

0.05 

2.91676 

18.00 

0.87108 

4 

20 

0.05 

2.91676 

19.00 

0.89107 

4 

20 

0.05 

2.91676 

20.00 

0.90835 

4 

20 

0.05 

2.91676 

21.00 

0.92321 

4 

20 

0.05 

2.91676 

22.00 

0.93591 

4 

20 

0.05 

2.91676 

23.00 

0.94672 

4 

20 

0.05 

2.91676 

24.00 

0.95586 

4 

20 

0.05 

2.91676 

25.00 

0.96356 
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NONCENTRALITY 

POWER 

NN  = 

5.0 

POWER  = 

0.30382 

NN  = 

6.0 

POWER  = 

0.36314 

NN  = 

7.0 

POWER  = 

0.42203 

NN  = 

8.0 

POWER  = 

0.47946 

NN  = 

9.0 

POWER  = 

0.53461 

NN  = 
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POWER  = 

0.58677 

NN  = 

11.0 

POWER  = 

0.63550 

NN  = 
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POWER  = 

0.68050 

NN  = 

13.0 

POWER  = 

0.72163 

NN  = 

14.0 

POWER  = 

0.75885 

NN  = 

15.0 

POWER  = 

0.79223 

NN  = 

16.0 

POWER  = 

0.82192 

NN  = 

17.0 

POWER  = 

0.84812 

NN  = 

18.0 

POWER  = 

0.87108 

NN  = 

19.0 

POWER  = 

0.89107 

NN  = 

20.0 

POWER  = 

0.90835 

NN  = 

21.0 

POWER  = 

0.92321 

NN  = 

22.0 

POWER  = 

0.93591 

NN  = 

23.0 

POWER  = 

0.94672 

NN  = 

24.0 

POWER  = 

0.95586 

NN  = 

25.0 

POWER  = 

0.96356 
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IV.  POWER  OF  THE  T-TEST 

A.  INTRODUCTION 

In  testing  the  hypothesis  that  the  mean  of  a  normal  distribution  is  equal  to  Hq 
when  the  standard  deviation  is  unknown,  the  t  statistic  may  be  used.  This  test-statistic 
is  a  function  of  the  sample  mean,  the  sample  standard  deviation,  and  the  sample  size. 

The  t-test  has  a  very  wide  application  in  research  areas.  For  example,  we  may 
wish  to  know  whether  product  A  is  better  than  product  B;  or  whether  the  outputs  of 
machines  C  and  D  form  a  homogeneous  mass  of  product;  and  so  forth.  To  answer 
such  questions  we  can  employ  the  t-test  which  is  one  of  the  most  useful  techniques  in 
the  field  of  statistical  inference. 

In  such  hypothesis-testing  situations  the  power  of  t-tests  is  usually  of  interest  to 
the  experimenter.  We  may  use  the  power  curve  to  evaluate  the  power  of  the  t-test. 
The  required  sample  size  might  be  determined  by  referring  to  a  set  of  power  curves. 
We  will  consider  the  power  curves  for  the  one-sample  and  two-sample  cases.   [Ref  10] 

B.  ONE-SAMPLE  T-TEST 

We  assume  that  a  random  sample,  XI,  X2,  X3,  ...,  Xn  of  size  n  is  taken  from  a 
normal  population  with  mean  ^  and  standard  deviation  <y.  Suppose  we  wish  to  test  the 
hypothesis  that  \i  <>  \Iq  against  the  alternative  hypothesis  that  fi  >  Hq  .  The 
procedure  is  to  reject  the  hypothesis  Ho  :  ji  <  ^q  if  ((X-  jIq  )  V  n  )  /  S  >  t  ^^  ,  and  to 
accept  Ho  if  otherwise  where  X  and  S  are  the  sample  mean  and  sample  standard 
deviation  and  where  t  ^  is  the  I  —  a  percentage  point  of  the  Student  t-distribution 
with  n-1  degrees  of  freedom,  i.e., 

P(  Student  t  >  t^  )  =  a  (4.1) 

With  critical  region  (1^,00),  the  hypothesis  ^  =  jIq  is  rejected  with  probability 
a  when  Ho  is  true.  If  the  mean  of  the  normal  distribution  is  fi^  >  JIq  ,  then  the 
power  of  the  above  test  is 

P(  noncentral  t  >  t  ^  |  6  =  ((  ^  j  -  fig)  Vn  )/  a)  =  1-  P  (4.2) 
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where  6  is  the  noncentrality  parameter  for  the  noncentral  t-distribution  with  n-1 
degrees  of  freedom  [Ref  6:p  .320]. 

C.  TWO-SAMPLE  T-TEST 

Suppose  we  are  given  two  random  samples  from  two  normal  populations  XI,  X2, 

...,  Xn  and  Yl,  Y2,  ...,  Ym,  where  the  X's  and  Y's  are  independent  with  mean  values 

7  7  0 

}i^  =  n    and  common  unknown  variances  <y^      =  <7        =  cr"-.    Suppose  we  wish  to 

test  the  null  hypothesis  Ho  :  n^  <  ^    against  the  alternative  hypothesis  Ha  that    ]i^ 

>  \iy.   The  procedure  is  to  reject  the  hypothesis  Ho  if  (X-  Y)  /  V  (  S^     (1/n  +  1/m)) 

>  t  ^  where  t  ^  is  a  critical  value  of  the  Student  t-distribution  based  on  m+n—  2 
degrees  of  freedom.  In  this  case  take  6  =  (  n^  —  fi  )  /  (T  V  (1/n  +  1/m)  so  the 
probability  of  rejecting  the  hypothesis  when  ^^  >  ^    is  equal  to 

P(  noncentral  t  >  t  ^  |  6  )  =   1-  p  (4.3) 

where  6  is  the  noncentrality  parameter  for  the  noncentral  t-distribution  with  m+n-2 
degrees  of  freedom  [Ref  6:p  .321]. 

D.  TWO-SIDED  T-TEST 

For  tests  where  the  alternative  hypothesis  specifies  that  yi  ¥=  n^  in  the  test  of 
section  B  or  ^^  t^  ji  .^  in  the  test  of  section  C,  where  the  direction  can  be  either  up  or 
down,  i.e,  ^  >  fi^  or  n  <  fig  and  \i^  >  \iy  or  Ji^  <  ^,.,  we  have  what  is  known  as  a 
two-sided  test.  In  this  case  the  rejection  procedure  for  the  test  of  section  B  specifies 
that  Ho  is  rejected  if  either  ((X-  jIq  )  V  n  )  /  S  >  t  q/2  or  ((X-  ji^  )  V  n  )  /  S  <  -  t  ^/2 
and  the  power  of  the  test  is 

P(  noncentral  t  >  t  qj/2  |  6  )  +  P(  noncentral  t  <  -t  ^/2  |  6  )  =   1  -  P  (4.4) 

where  6  =  (  fij  -jIq  )  V  n  /  (T  [Ref  6:p  .323]. 


E.       THE  ALGORITHMS  AND  FLOWCHART 

Listings  of  programs  for  generating  the  power  table  and  curve  for  t-tests  are 
given  in  the  Appendix  E. 

The  following  steps  are  used  to  compute  the  power: 
1)     Determining  whether  one-sample  or  two-sample  t-tests  should  be  used. 
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2)  After  determining  the  test,  give  the  input  data,  (n,  a  ,  ^g  ,  a  ) 

3)  Determine  whether  the  one-sided  or  two-sided  test  should  be  used. 

4)  Compute  the  critical  region  (inverse  t-values).    (The  computation  method  is 
provided  in  Appendix  B) 

5)  Compute    the    power    corresponding    to    the    chosen    user    option.      (The 
computation  method  is  discussed  in  Chapter  II) 

6)  Plot  the  power  corresponding  to  each  above  option. 
A  flowchart  is  shown  in  Figure  4.1a 

F.       EXAMPLES  OF  T-TESTS 

1.  Example  of  one  sample  one-sided  t-test 

a.  Scenario 

The  supervisor  in  an  electronic  instrument  plant  is  concerned  with  the 
repair  time  of  an  electronic  equipment.  Standards  specify  that  the  mean  repair  time 
must  be  at  most  40  hours.  For  power  computation,  it  is  assumed  that  the  variance  of 
repairing  time  is  4.0.   The  appropriate  hypothesis  are 

Ho  :  ^  >  40, 

Ha  :  n  <  40. 
The  experimenter  decides  to  use  a  sample  of  n=  25  observations  and  a  =  .05.  He  would 
like  to  know  what  the  power  of  the  t-test  is  for  a  range  of  possible  population  means. 
He  decides  to  check  means  from  38  to  40  in  .1  increments. 

b.  Inputs 


INPUTS:  step  1 

step  2 

step  3 

step  4 

step  5 

step  6 
c.   Start  program 
START  PROGRAM: 
Do  you  want  to  test  one  sample  t-test(y/n)? 
y 
Sample  size  ? 

25 
Alpha-level  ? 


one  sample  one-sided  t-test 

Mo  =  40 

(T  =  2.0 

a  =.05 

sample  size  =  25 

maximum  Hj  =40.,  minimum  Hj  =38,  increment  =  .1 
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y 
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jne 


-VL 


Ho   :  I-fJ  =  rfJzero 
Ha   :   m  >   MLIzero 


V 


$ 


""^      )j         Ta' 


/     ,- 


o  ^  i  ci'^  '^  -^  ^  -  —   0 


i_ 


Ho   :  WJ  =  I-rJzero       /    !J    . 


Ha   :   MU  <  r-fJzero 


c 


i_ 


Stop 


■7 
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Figure  4.  la     System  flowchart  for  t-test  power. 
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Figure  4.1b  System  nowchart  for  t-iest  power(cont'd.) 
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Figure  4.1c  System  flowchart  for  t-test  power(cont'd.) 
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Figure  4. Id  System  flowchart  for  t-test  power(cont'd.) 
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? 

.05 
Mu-zero  value  ? 

40 

Standard  deviation  value  ? 

7 

2 

Do  you  want  to  use  one-sided  test(y/n)? 

y 

Do  you  want  to  test  Ho  :  mu    ^  muO  vs  Ha  :  mu  >  muO  (y/n)? 

n 

Do  you  want  to  test  Ho  :  mu    ^  muO  vs  Ha  :  mu  <  muO  (y/n)? 

y 

Maximum  mul  value  (the  maximum  value  on  X-axis)? 

? 

40 

Minimum  mul  value  (the  minimum  value  on  X-axis)? 

(condition  :  mul  must  be  less  than  muO  ) 

7 

38 

Increment  mul  value  ? 

7 

.1 

d.   Output 

The  screen  output  is  as  follows: 
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h 

] 

30wer 

NN  = 

38.000 

POWER  = 

0.99933 

NN  = 

38.100 

POWER  = 

0.99849 

NN  = 

38.200 

POWER  = 

0.99677 

NN  = 

38.300 

POWER  = 

0.99346 

NN  = 

38.400 

POWER  = 

0.98742 

NN  = 

38.500 

POWER  = 

0.97706 

NN  = 

38.600 

POWER  = 

0.96028 

NN  = 

38.700 

POWER  = 

0.93465 

NN  = 

38.800 

POWER  = 

0.89772 

NN  = 

38.900 

POWER  = 

0.84755 

NN  = 

39.000 

POWER  = 

0.78326 

NN  = 

39.100 

POWER  = 

0.70559 

NN  = 

39.200 

POWER  = 

0.61711 

NN  = 

39.300 

POWER  = 

0.52203 

NN  = 

39.400 

POWER  = 

0.42564 

NN  = 

39.500 

POWER  = 

0.33355 

NN  = 

39.600 

POWER  = 

0.25053 

NN  = 

39.700 

POWER  = 

0.17996 

NN  = 

39.800 

POWER  = 

0.12338 

NN  = 

39.900 

POWER  = 

0.08062 

40 


0.99933 
power 


^_l<t__)<t J'C  —  —  *^.  —  —  *  —  —  *  —  —  —  ^—  —  . 


■+-+ 


0.53997 


0.  80615E-01     +-+ + + + 

38.00  38.95 

X-SCALE:   "-"=  0.237E-01  UNITS 


\i,         39.90 


Y-SCALE: 


=  0.153E-01  UNITS 


2.  Example  of  One  sample  two-sided  t-test 
a.   Scenario 

The  supervisor  in  a  rocket  propellant  plant  is  concerned  with  the  burning 
rate  of  a  rocket  propellent.  Standards  specify  that  the  mean  burning  rate  must  be  40 
inches  per  second.  For  power  computation,  it  is  assumed  that  (T  =  4.0.  The 
appropriate  hypotheses  are 

Ho  :  fi  =  40 

Ha  :  ^  #  40 
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INPUTS: 


The  engineer  decides  to  use  a  sample  of  n=  25  observations  and  a  =  .05.  He  would  like 
to  know  the  power  of  the  t-test  is  for  a  range  of  possible  population  means.    He 
decides  to  check  means  from  38  to  42  in  .1  increments. 
b.   Inputs 

one  sample  two-sided  t-test 

Ho  =  40 

<T  =  2.0 

a  =.05 

sample  size  =  25 

maximum  flj  =42,  minimum  Hj  =38,  increment  =.1 


step  1: 

step  2: 

step  3: 

step  4: 

step  5: 

step  6: 
c.   Start  program 
START  PROGRAM: 
Do  you  want  to  test  one  sample  t-test(y/n)? 

y 

Sample  size  ? 
7 

25 
Alpha-level  ? 

7 

.05 
Mu-zero  value  ? 

7 

40 

Standard  deviation  value  ? 

7 

2 

Do  you  want  to  use  one-sided  test(y/n)? 

n 

Do  you  want  to  use  two-sided  test(y/n)? 

y  -**•■■ 

Maximum  mul  value  (the  maximum  value  on  X-axis)? 

7 

42 

Minimum  mul  value  (the  minimum  value  on  X-axis)? 

(condition  :  mul  must  be  less  than  mu-zero) 
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7 

38 

Increment  mul  value  ? 

7 

.1 

d.    Output 

The  screen  output  is  as  follows: 


POWER 


NN  = 

38.000 

POWER  = 

0.99770 

NN  = 

38.100 

POWER  = 

0.99525 

NN  = 

38.200 

POWER  = 

0.99073 

NN  = 

38.300 

POWER  = 

0.98279 

NN  = 

38.400 

POWER  = 

0.96965 

NN  = 

38.500 

POWER  = 

0.94910 

NN  = 

38.600 

POWER  = 

0.91875 

NN  = 

38.700 

POWER  = 

0.87639 

NN  = 

38.800 

POWER  = 

0.82057 

NN  = 

38.900 

POWER  = 

0.75109 

NN  = 

39.000 

POWER  = 

0.66944 

NN  = 

39.100 

POWER  = 

0.57882 

NN  = 

39.200 

POWER  = 

0.48379 

NN  = 

39.300 

POWER  = 

0.38975 

NN  = 

39.400 

POWER  = 

0.30192 

NN  = 

39.500 

POWER  = 

0.22459 

NN  = 

39.600 

POWER  = 

0.16066 

NN  = 

39.700 

POWER  = 

0.11146 

NN  = 

39.800 

POWER  = 

0.07708 

NN  = 

39.900 

POWER  = 

0.05692 

NN  = 

40.000 

POWER  = 

0.05028 

NN  = 

40.100 

POWER  = 

0.05687 

NN  = 

40.200 

POWER  = 

0.07698 

NN  = 

40.300 

POWER  = 

0.11131 

NN  = 

40.400 

POWER  = 

0.16045 
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NN  = 

40.500 

POWER  = 

0.22433 

NN  = 

40.600 

POWER  = 

0.30162 

NN  = 

40.700 

POWER  = 

0.38942 

NN  = 

40.800 

POWER  = 

0.48344 

NN  = 

40.900 

POWER  = 

0.57848 

NN  = 

41.000 

POWER  = 

0.66912 

NN  = 

41.100 

POWER  = 

0.75081 

NN  = 

41.200 

POWER  = 

0.82034 

NN  = 

41.300 

POWER  = 

0.87621 

NN  = 

41.400 

POWER  = 

0.91862 

NN  = 

41.500 

POWER  = 

0.94901 

NN  = 

41.600 

POWERS 

0.96959 

NN  = 

41.700 

POWER  = 

0.98275 

NN  = 

41.800 

POWER  = 

0.99070 

NN  = 

41.900 

POWER  = 

0.99524 
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V.  DESCRIPTION  OF  THE  POWER  PROGRAM 

A.  PROCEDURE  OVERVIEW 

The  interactive  program  included  in  this  thesis  is  written  in  FORTRAN  77.  It 
was  written  for  u:e  on  an  IBM370  from  IBM3278  terminal.  It  is  an  interactive 
program.  The  command  'NONCENT'  is  all  that  is  required  to  start  this  program.  The 
executive  file  name  initializes  the  virtual  machine  environment  and  asks  the  user 
questions  about  the  choice  of  program  and  compilation  requirements.  It  then  activates 
the  selected  program.  Both  ANOVA  and  TTEST  are  interactive  programs  which  will 
be  discussed  in  detail  in  later  chapters.  The  output  from  the  programs  is  presented  on 
the  terminal  screen. 

B.  INSTRUCTIONS  FOR  PROGRAM  ACCESS 

To  start  the  program,  make  sure  you  have  loaded  (NONCEN,  ANOVA,  and 
TTEST)  on  your  disk.  In  CMS  (operating  system  mode),  type  'NONCEN'  and  you  will 
see: 

Please  provide  the  FILENAME  for  your  VS  FORTRAN  program. 
Now  type  the  program  name  you  want,  for  example  'TTEST',  the  response  is: 

Do  you  need  to  compile  your  program  ?  (y/n) 
If  you  want  to  run,  type  'Y'  and  your  program  will  be  loaded. 

A  detailed  description  of  the  ANOVA  (Analysis  of  Variance  Test  )  is  given  in 
Chapter  III  and  of  the  t-test  in  Chapter  IV. 

Following  the  screen  output  from  the  selected  test,  the  user  will  be  asked  some 
questions  about  the  output: 
Do  you  wish  to  BROWSE  your  output?  (Y) 
n 

Print  your  output  file?  (Y) 
n 

Do  you  wish  to  XEDIT  the  program  file?  (Y/N) 
n 

Do  you  wish  to  run  the  program  again?  (Y) 
n 

Then  return  to  CMS  mode. 
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VI.  SUMMARY  AND  CONCLUSIONS 

Power  considerations  are  useful  in  the  design  and  assessment  of  statistical  tests. 
The  calculation  of  power  usually  involves  noncentral  distributions  for  which  tables  of 
probability  are  not  available.  A  search  was  made  for  algorithms  approximating  the 
noncentral  t,  F  and  x  distributions.  Algorithms  giving  sufficient  accuracy  and 
making  efficient  use  of  computer  resources  have  been  implemented  in  this  thesis. 
Listings  of  the  FORTRAN  code  for  these  implementations  are  included. 

An  interactive  program  to  compute  and  display  power  curves  for  several  t-test 
and  F-test  situations  has  been  developed.  This  program  is  user  friendly  and  is 
described  in  this  thesis.  A  listing  of  the  program  is  provided.  This  program  should  be 
useful  to  researchers,  experiment  designers  and  statisticians. 
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APPENDIX  A 
AN  INVERSE  CENTRAL  F  APPROXIMATION 

Let     P     be     the     central     F    complementary    CDF,     and     let    Q     be    the 
complementary  standard  normal   CDF  [Ref.    4:  p.    947]. 
Define  yp  by  Q(yp)=p,   and  Fp  by  P(Fp   |    v^   ,V2  )=p. 
Then       Fp   ~   exp(2w) 
where     w  =  wl-(w2)(w3) 

wl  =  (yp  V  (h+X))/h 

w2  =  (1/(V2  -l)-l/(vi  -1)) 

w3  =  (X+(5/6)+2/(3h)) 

h  =  2/((l/(vi  -1)+1/(V2  -1)) 


2  .- 
The  program  we  used  in  the  computation  of  critical  values  of  F-tests  is  shown 


X=(yp  ^  -3)/6 


below. 

V  u  ^^^  t^\f\t\^%^\^\ftj  ^^  tj  \j  \j  \j  s^  ^f  u  tj  \j  •t*\^'v\£s^\f\y\jsf\^\^\^  "i^  \j  ^^  %^  \f  %t  \r  t^  ^  \f  \^  \f  \f  \^  \£  \*  \y  \f  ^^\j\f\f\t\^\j\f\f\^  xr  \rf  x^  ■^  vf  \^  vf 

^  W  TWTT  W  K  K  K  H  K  K  KW  ^  W  K  K  XK^K  W  WW  WWWWAWRA^AAR  W  K  K  K  W  W  W  R  W  K  W  W  K  W  W  W  W  W  XXMAXXXXXK  W  W  K  K  X  R  K 

N  W 

*  APPROXIMATION  TO  THE  F-  INVERSE  DISTRIBUTION.  * 

*  CALCULATE  F-INVERSE  GIVEN  DOFl ,D0F2 , ALPHA.  » 
»       Q(  F-ALPHA  |D0F1,D0F2  )=  ALPHA.  » 

*  K 

KWWWWWW  AW*  .  .  *W*W:ffWW**WKWWWW*WWWW"^*WWWWWWW!W*WWWW9T*WWWWWKWWWWWWWWWWWW*WWWW 

REAL  ALPHA,NUMER,DENUM,LAM,XP,YP,FINVER,A,B 

DATA  NF1,NF2, ALPHA/2, ^,.01/ 

DATA  CO ,C1 ,C2/2 . 515517 ,0 . 802853 ,0 . 010528/ 

DATA  Dl, 02, D3/l.«2788, 0.189269, 0.001508/ 

B=NFl/2. 

A=NF2/2. 
C 

C     CALCULATE  NORMAL  QUANTILE. 
C 

T=SQRT1AL0G(1./(ALPHA*»2. ))) 

NUMER=C0+C1*T+C2*T**2 

DENUM=1 . +D1*T+02»T»*2 . +D3*T**3 . 

XP=T-( NUMER/DENUM ) 
C 

C     CALCULATE  THE  INVERSE  F-DISTRIBUTION. 
C 

LAM=(XP»«2-3. )/6. 

M2=  (l./(2.*B-l.  ))-(!./( 2. »A-1. )) 

W<i=  (l./(2.»A-l.  ))  +  (!./(  2. *B-1.  )) 

H=2.*m««(-1.  ) 

Wl  =  (LAM*(5./6.  )-(2./(3.*H))) 

W5=(XP»(H+LAM)**.S)/H 

W=W3-H2»W1 

FINVER=EXP(2»W) 

WRITE(5,100)  NF1,NF2, ALPHA, FINVER 
100   FORMAT( '0','   NF1=' ,I5,'NF2=' ,15 , * ALPHA= ' ,F5.4, ' FINVERSE= ' ,F6.3) 

STOP 

END 

The  following  table  shows  some  values  obtained  using  this  approximation, 
together  with  exact  values  of  critical  values  of  the  F  distribution.  It  can  be  seen  in 
Table  8  that  this  approximation  provides  adequate  accuracy  for  the  computation  of  the 
critical  value  of  the  F-tests. 
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TABLC  8 
APPROXIMATION  AND  HXACT  CRITICAL  VALUHS  OF  F-  ins  IS 


a 

nl 

n2 

approx 

exact 

.05 

4 

4 

6.581 

6.39 

4 

6 

4.60 1 

4.53 

4 

10 

3.525 

3.48 

10 

20 

2.35 

2.35 

20 

30 

1.932 

1.93 

6 

4 

6.359 

6.16 

10 

6 

4.112 

4.06 

20 

6 

3.944 

3.87 

30 

20 

2.041 

2.04 

30 

30 

1.841 

1.84 

.10 

4 

4 

4.150 

4.11 

4 

6 

3.218 

3.18 

4 

10 

2.648 

2.61 

10 

20 

1.939 

1.94 

20 

30 

1.668 

1.67 

6 

4 

4.046 

4.01 

10 

6 

2.952 

2.94 

20 

6 

2.864 

2.84 

30 

20 

1.739 

1.74 

30 

30 

1.607 

1.61 
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APPENDIX  B 
THE  INVERSE  CENTRAL  T  APPROXIMATION 

We  can  use  the  asymptotic  expansion  for  the  Inverse  t-CDF,  as 
follows.  Let  Q  be  the  complementary  CDF  of  standard  normal,  and  let  P 
be  the  central  t  CDF  [Ref.  4:  p.  949]. 


then 


If  P(tp  I  V  )  =  l-2p  and  Q(Xp  )  =  p 


tp^  Xp  +  (91  (Xp  )/  v)+  (92  (Xp  )  /  v2  )-H(g3  (Xp  )  /  v^  )+(  g^  (Xp  )  / 


where  9-^  (x)=.  25(x-^  +x), 

92  (x)=(l/96)(5x^  +16x^  +3x), 

93  (x)=(l/384)(3x^  +19x^  +17x^  -15x),  and 

94  (x)=(l/92160)(79x^  +776x^  +1482x^  -1920x^  -945x). 
The  program  we  used  in  the  computation  of  the  critical  value  of  the 
t-test  is  shown  below. 

« 

»       APPROXIMATION  TO  THE  T-  INVERSE  DISTRIBUTION 
»       CALCULATE  T-INVERSE  GIVEN  DOF  ALPHA. 
*       Q(  T-ALPHA  I  DOF  )=  ALPHA. 
* 

WWWW^WKWWARWWWWAAAA^W^^AAA^WAAWXAA^AAWAKAAAAAAAKAAARAKAAAAA 

c 

C     set  initial  conditions. 
C 

REAL  ALPHA ,NUMER ,DENUM > LAM ,XP >G1 ,62 ,GZ ,X,TINVER 

INTEGER  N 

DATA  N, ALPHA/9,. 025/ 

DATA  CO, CI, 02/2.515517,0. 802853, 0.010228/ 

DATA  D1,D2, 03/1.432788,0. 189269, 0.001508/ 


100 


CALCULATE  NORMAL  QUANTILE. 

ALPHA=.05 

T=SQRT(AL0G(1./(ALPHA**2.  ))) 

NUMER=CC+C1*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP=T-(NUMER/DENUM) 

CALCULATE  THE  INVERSE  T-DISTRIBUTION. 

X=XP 

G1=.2S*(X**3+X) 

G2  =  (l . /96 .  )*( 5 . *X**5+16 . *X**3+3 . »X ) 

G3  =  (l./38<*.  )*(3.»X**7+19.*X**3+17.»X»*3-15.*X) 

G4=( 1./92160.  )*( 79.*X**9+776.*X**7*1<+82.*X**5-1592.»X**3-945.»X) 

TP=X+(G1/N)  +  (G2/(N*N))  +  (G3/(N»*3))  +  (G'+/(N**4)) 

TINVER=TP 

HRITE(S,100)  N, ALPHA, TINVER 

FORMAK  'OS'   D0F=',I3,4X,'ALPHA=',F5.3,6X,'TINVERSE=',F9.5) 

STOP 

END 


The  following  table  shows  some  values  obtained  using  this 
approximation  together  with  exact  values  of  critical  values  of  the  t 
distribution.   It  can  be  seen  in  Table  9  that  this  approximation 
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provides  adequate  accuracy  for  computing  the  critical  regions  of  the 
t-tests. 

TABLE  9 
APPROXIMATION  AND  EXACT  CRITICAL  VALUES  OF  T-TESTS 


DOF 

a  (  Type 

I  error  ) 

.  05 

.  01 

exact 

approx 

exact 

approx 

2 

2.  919 

2.  867 

6.  964 

6.  475 

4 

2.  132 

2.  126 

3.  746 

3.  701 

6 

1.  943 

1.  942 

3.  142 

3.  130 

8 

1.  833 

1.  833 

2.  821 

2.  818 

10 

1.  812 

1.  823 

2.  764 

2.  761 

12 

1.  783 

1.  783 

2.  681 

2.  680 

14 

1.  761 

1.  761 

2.  624 

2.  624 

16 

1.  746 

1.  746 

2.  583 

2.  583 

18 

1.  734 

1.  734 

2.  552 

2.  552 

20 

1.  725 

1.  725 

2.  528 

2.  528 

25 

1.  708 

1.  708 

2.  485 

2.  485 

35 

1.  690 

1.  690 

2.  438 

2.  438 

45 

1.  690 

1.  690 

2.  412 

2.  413 

65 

1.  669 

1.  669 

2.  385 

2.  385 

85 

1.  663 

1.  663 

2.  371 

2.  371 
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APPENDIX  C 
AN  APPROXIMATION  OF  THE  NORMAL  CDF 

Let  Z  be  the  standard  normal  density  function  and  let  P  be  the 
corresponding  CDF  [Ref.  4:  p.  939].   Then 

P(X)  =  J  Z(t)dt  -   l-Z(x)(ait+a2t2+a3t3)+  e  (x) 

where  t=l/(l+px) 

I  £  (x)|  <  .000001 

p  =.33267,  and         a^  =  .4361836 

32  =  -  . 1201676 

33  =  .9372980 
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APPENDIX  D 
THE  COMPUTATION  OF  NONCENTRAL  F  DISTRIBUTION 


»^»»»»»*jf»»)0(KXXXXXXKXj^   INFORMATION   ^^^^»^^^^»»»»^H^^»»y  ^»»^^x^»»  x  x  x  x  »» 

»        THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  F    * 

*  DISTRIBUTION  ( CDF,F-INVERSE  ) .  * 
»  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  » 
w  * 
*»^HHH«f****a*********»*  VARIABLE  DEFINITION  *************************** 

*  * 

*  FVAL  :  X-VALUE  DIMENSION  » 
»  FCDF  :  NON-CENTRAL  F  CDF  DIMENSION.  * 

*  FINVER  :  GIVEN  X  VALUE,  FIND  INVERSE  VALUE.  * 
»  FLAM  :  NON-CENTRAL  PARAMETER  * 
»  FX  :  TO  FIND  INVERSE,  GIVEN  X  VALUE.  » 

*  FSTAR  :  X  VALUE  VARIABLE.  * 
»  NFI  :  DEGREE  OF  FREEDOM  1.  * 
«  NF2  :  DEGREE  OF  FREEDOM  2.  * 

*  CDF  :  NORMAL  DISTRIBUTION  FUNCTION.  * 
»  P  (  Z  <=  X)  * 
«  CONST  :  1  /  SORT  (  2  *  PI  )  IN  NORMAL  DISTRIBUTION  FUNCTION.  * 
«  B  :  THE  VALUE  OF  'X'  IN  NORMAL  DISTRIBUTION.  * 
«  * 

A  A  K  A  A  A!  A  A  AAA  A  AAA  AAA  W*  W  AT  A  ^W  WW  AAAAAAAAA  W^  A^  W  A  A  A  A  A  A  W  WWW  W  A  A  A  A  WWW  WW  AAAAAAAAAAA 

C 

C     SET  INITIAL  CONDITIONS. 
C 

REAL  FVAL(100),FCDF(100), FINVER, FLAM, FMAX,FX,FINIAL 
DATA  NFl,NF2,FLAM,FMAX/2,4,75.0,90./ 
DATA  Al  ,A2  ,A5 ,P/.«61836 ,- .  1201676 , .  9372980 , .  32267/ 
DATA  FSTAR/25.128/ 

WRITE(3,2<i)  NFI, NF2, FLAM 
WRITE! 3,21) 
WRITE(3,22) 
WRITE! 3,23) 

CONST   =  1.0  /  SQRT  (2.0  *  3.1415927) 
C 

C     CALCULATE  THE  NORMAL  APPROXIMATION  . 
C 

Cl=( NF1*FSTAR  )/( NFl+FLAM ) 

C2=l-(2./(9.*NF2)) 

C3=(2.*(NF1+2.*FLAM))/(9.*(NF1+FLAM)»*2. ) 

CCl=(Cl»*(l./3.  ))»C2-(1-C3) 

CC2=((C3+(2./{9.»NF2))*(Cl)**(2./3. ))**.5) 

FX=CC1/CC2 

CDF  =  0.0 

1F(FX  .EQ.O)  THEN 

CDF  =  0.5 
ELSE  IF(FX.GT.O)  THEN 
T=1./(1.+P*FX) 
XSQU=(FX*FX)*.S 
2PDF=C0NST*EXP( -XSQU ) 
CDF=1-ZPDF*( A1*T+A2*T»*2+A3*T**3 ) 
ELSE 
FX=-FX 

T=1./(1.+P*FX) 
XSQU=(FX»FX)«.5 
ZPDF=CONST*EXP( -XSQU) 
CDF=ZPDF*( A1*T+A2*T**2+A3»T»*3 ) 
END  IF 

FVAL ( FSTAR )= FSTAR 
FCDF( FSTAR )=CDF 
HRITEl 3,250)  FSTAR, , FCDF ( FSTAR ) 
250      FORMAK  '0',F7.3,8X,F8.5) 

21  FORMATC  •0',«*0(  '_' )) 

22  FORMAT( ■0',2X,"'X"-VALUE',8X,'F(Z  <=  X)') 

23  FORMAT!  'OS^OC  '_•  )) 

2<f       FORMAT! '0','DOF1=', 12, 3X,'DCF2=', 12, 3X, 'NON-CENTRAL  PARAMETER= 
8■,F8.c^) 
STOP 
END 
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APPENDIX  E 
THE  COMPUTATION  OF  NONCENTRAL  T  DISTRIBUTION 


*********»*»»******»*»»   INFORMATION   ^KJtKXXXXKJHtXKKKXKMXXMKKXXXKKXXXKXX 

*  * 

»  THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  T  * 

*  DISTRIBUTIONS  (  CDFs )  * 

*  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  * 

*  (REF.  SEVERO  8  ZELEN'S  (1960)  )  » 

*  * 
**********************  VARIABLE  DEFINITION  »**»xxxMXitMXJtxxxxxxxxKX*x*» 

*  ^ 

*  TVAL  :  X-VALUE  DIMENSION  * 

*  TCDF  :  NON-CENTRAL  T  CDF  DIMENSION.  » 

*  TLAM  :  NON-CENTRALITY  PARAMETER  » 

*  TX  :  TO  FIND  INVERSE,  GIVEN  X  VALUE.  * 

*  TSTAR  :  X  VALUE  VARIABLE.  * 

*  NF   :  DEGREE  OF  FREEDOM  * 

*  CDF   :  NORMAL  DISTRIBUTION  FUNCTION.  * 

*  P  (  Z  <=  X)  * 

*  CONST  :  1  /  SORT  (  2  *  PI  )  IN  NORMAL  DISTRIBUTION  FUNCTION.   * 

*  B      :  THE  VALUE  OF  'X'  IN  NORMAL  DISTRIBUTION.  * 

c 

c     SET  INITIAL  CONDITIONS. 

c 

REAL  TVAL(  100 ) ,TCDF( 100 ) ,TINVER,TLAM,TMAX,TX,TINIAL 
DATA  NF,TLAM/24,53.078/ 

DATA  A1,A2, A3, P/. '^361836,-. 1201676,. 9372980,. 33267/ 
DATA  TSTAR/24.<i9^9/ 
WRITE (3,2'^)  NF,TLAM 
WRITE(3,21) 
WRITE(3,22) 
WRITE(3,23) 

CONST   =  1.0  /  SORT  (2.0  »  3.1^15927) 
c 

C     CALCULATE  THE  NORMAL  APPROXIMATION  . 
c 

Cl=( 1. -( l./(4.*NF ) ) )*TSTAR-TLAM 
C2=SQRT( l.+(TSTAR**2/( 2.*NF  ) ) ) 
TX=C1/C2 

CDF  =  0.0 
IF(TX  .EQ.O)  THEN 
CDF  =  O.S 

ELSE  IF(TX.GT.O)  THEN 
T=1./(1.+P*TX) 
XSQU=(TX»TX)*.S 
ZPDF=CONST*EXP( -XSQU ) 
CDF=1-ZPDF*( A1*T+A2*T»*2+A3«T*»3 ) 
ELSE 
TX=-TX 

T=1./(1.+P*TX) 
XSQU=(TX*TX)».5 
ZPDF=CONST»EXP( -XSQU ) 
CDF=ZPOF*( A1*T+A2»T**2+A3»T»*3 ) 
END  IF 

TVAL ( TSTAR )=TSTAR 
TCDF( TSTAR )=CDF 
WRITE (3,250)  TSTAR, , TCDF (TSTAR ) 
'0',F7.3,8X,F8.5) 
■0',40(  •_'  )) 

'0',2X,"'X"-VALUE',8X,'F(Z  <=  X)') 
'0','»0(  ■_• )) 
'0','DOF=' , I 2, 3X, 'NON-CENTRAL  PARAMETER= 


250 

FORMAT ( 

21 

FORMAT ( 

22 

FORMAT ( 

23 

FORMAT ( 

24 

FORMAT ( 

a- 

,F8.4) 

STOP 

END 
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APPENDIX  F 
THE  POWER  OF  ANOVA  PROGRAM  LIST 


POHER  ALGORITHM  (  ANALYSIS  OF  VARIANCE  ) 


DIRECTED  BY  :  PROFESSOR  DONALD. R.BARR 


JULY  1986 
NAVAL  POSTGRADUATE  SCHOOL 


* 

********************'***Jf******^t*»****»*^Hf***if  ************************** 

REAL  NN{ 200  ),POWER( 200  ) , NON,FINV,NONCEN,FPOW,FX, ALPHA, CDF , CONST 

REAL  MAXDEL ,MINDEL ,INCDEL ,STDS, PICON, MINSD ,MAXSD ,INCSD ,MAXALP 

REAL  MINALP,INCRAL 

INTEGER  MF( 20  ),MAXN, MINN, INCRN,MAXK, MINK, INCRK,JJ,K,N 

CHARACTER*!  ANS 

C 

C 
C 
C 

C 

500 


* 
* 
* 
* 
* 
* 

*  WRITTEN  BY  :  HUR,  SEONG  PIL 

*  DEPARTMENT  OF  O.R 
* 

* 

*  PROGRAM  USES  THE  POWER  PLOT  TO  TEST  FOR  ANALYSIS  OF  VARIANCE 
* 

* 

* 
* 

* 

* 
* 

*  NOTE 
* 


(A)  ONE-WAY  ANOVA 

(1)  SAMPLE  SIZE  VS  POWER 

(2)  NUMBER  OF  TREATMENT  VS  POWER 

(3)  ALPHA-LEVEL  VS  POWER 

(<♦)   NONCENTRALITY  PARAMETER  VS  POWER 

(B)  MULTI-WAY  ANOVA 

( 1  )   NUMBER  OF  REPLICATES  VS  POWER 

(3)  ALPHA-LEVEL  VS  POWER 

(4)  NONCENTRALITY  PARAMETER  VS  POWER 
subroutine  PLOTT  is  NONIMSL  subroutine  library. 


C_ 

C 

C 

C 

C 

C_ 

10 


START  ONE-WAY  OR  MULTI-WAY  ANALYSIS  OF  VARIANCE, 


PRINT* , ■ DO  YOU  WANT 

READ(5,7)  ANS 

FORMAT(Al) 

IF( (ANS.EQ. 'N' ).0R.( 

PRINT*, 'DO  YOU  WANT 
«  VS  POWER  (  Y/N)?' 

READ(5,7)    ANS 

IFKANS.EQ.  'Y'  ).0R.( 

PRINT*, 'DO  YOU  WANT 
S    (    Y/N)?' 

READ(5,7)    ANS 

IFdANS.EQ.  "Y"  ).0R.( 

PRINT*, 'DO  YOU  WANT 

READ(S,7)    ANS 

IF((ANS.EQ. 'Y' ).0R.( 

PRINT*, 'DO  YOU  WANT 

READ(5,7)    ANS 

IF((ANS.EQ. 'Y' ).0R.( 

PRINT*, 'DO  YOU  WANT 

READ(5>7)   ANS 

IFKANS.EQ.  'Y*  ).0R.( 

GO  TO  999 


TO  USE  ONE-WAY  ANOVA  (Y/N)?' 


ANS.EQ. 'N'  ))  GO  TO  1000 

TO  PLOT  N(  n   OF  OBSERVATION  PER  TREATMENT) 


ANS.EQ. 'Y' ))  GO  TO  10 
TO  PLOT  K  (  NUmER  OF 


TREATMENT  )  VS  POWER 


ANS.EQ. 'Y' ))  GO  TO  20 

TO  PLOT  ALPHA-LEVEL  VS  POWER  (Y/N)?' 

ANS.EQ. 'Y' ))  GO  TO  30 

TO  PLOT  NONCENTRALITY  VALUE  VS  POWER  (Y/N)?' 

ANS.EQ. 'Y' ))  GO  TO  40 
TO  PLOT  AGAIN  (Y/N)?' 

ANS.EQ. 'Y' ))  GO  TO  500 


ANALYZE  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  ONE-WAY 
ANALYSIS  OF  VARIANCE. 


PRINT*, 'NUMBER  OF  K(  »  OF  TREATMENT  )?' 

READ  (5,*)  K 

PRINT*, 'ALPHA-  LEVEL?' 

READ  (5,*)  ALPHA 

PRINT*,'  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  (5,*)  STDS 
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COMPUTE 

CENTRAL 

F-INVERSE  GIVEN 

D0F1,D0F2 

AND 

ALPHA-LEVEL. 

FINV=FINVER  (NFl 

,NF2„ 

ALPHA) 

COMPUTE 

NONCENTRALITY 

PARAMETER 

AS 

SAMPLE 

SIZE 

CHANGE . 

NONCEN=l 

^*STDS 

PRINT*,'  MAXIMUM  N  VALUE  (  THE  MAXIMUM  VALUE  ON  X-AXIS  )?' 

READ  (5,»)  MAXN 

PRINT*, 'MINIMUM  N  VALUE  (  THE  MINIMUM  VALUE  ON  X-AXIS  )? 
8  (CONDITION  :N  MUST  BE  MORE  THAN  2  )' 

READ  (5,»)  MINN 

PRINT*, 'INCREMENT   N  VALUE  ?' 

READ  (5,*)  INCRN 

WRITE! 6,95) 
95    FORMAT! IX, 'DOFl   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POHER') 

NF1=K-1 

JJ  =  1 

DO  110  N=MINN, MAXN, INCRN 

NF2=K*(N-1) 

C 

C 
C 
C 
C 

C 
C 
C 
C 

C 

c_ 
c 

C      COMPUTE  POWER  USING  NONCENTRAL  F-DISTRIBUTION  (CDF)  GIVEN  DOFl, 

C        D0F2, CENTRAL  F-INVERSE  AND  NONCENTRALITY  VALUE. 

C 

c 

POW=FPO( ALPHA, NFl, NF2,FINV,N0NCEN) 

WRITE ( 6 , 15  )  NFl ,NF2 , ALPHA , FINV ,NONCEN ,POW 
15    FORMAT!  '0 ' ,I5,SX,I5,2X,2F10 .S,1X,F12 . 2,1X,F10.5  ) 

NN(  JJ)=N 

POWER! JJ)=POW 

JJ=JJ+1 
110   CONTINUE 

JJ=JJ-1 

PRINT*, 'NN        POWER' 

DO  76  N=1,JJ 

WRITE(6,75)   NN( N), POWER! N) 

75  FORMAT! '0',  'NN= ' ,F8.3,^X, 'POWER= ' ,F10.5 ) 

76  CONTINUE 

C 

C 

C     PLOT  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  ONE-WAY 

C       ANALYSIS  OF  VARIANCE  USING  LIBRARY  SUBROUTINE. 

C 

c 

CALL  PLOTT!NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  !Y/N)?' 

READ  !5,7)  ANS 

IF!!ANS.EQ. 'Y' ).OR.!ANS.EQ. 'Y' ))  GO  TO  500 

GO  J?.    999 

C 

C 

C   ANALYZE  NUMBER  OF  TREATMENT  VS  POWER  IN  ONE-WAY  ANALYSIS  OF  VARIANCE 

C 

C 

20    PRINT*, 'tt  OF  OBSERVATIONS  PER  TREATMENT  ! N=  )  ?' 

READ  !5,*)  N 

PRINT*, 'ALPHA  LEVEL  (REAL  VALUE)  ?' 

READ  !5,*)  ALPHA 

PRINT*, 'WHAT  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  !S,*)  STDS 

PRINT*, 'MAXimJM  K  !»  OF  TREATMENTS)  RANGE  ?' 

READ  (5>N)  MAXK 

PRINT*, 'MINIMUM  K  (»  OF  TREATMENTS)  RANGE  (K  MUST  BE  MORE  THAN  2)?' 

READ  (5>*)  MINK 

PRINT*, 'INCREMENT   K  VALUE  ?' 

READ  (5,«)  INCRK 

WRITE(6,19S) 
195    F0RMAT(1X,'D0F1   D0F2  ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER') 

JJ=1 

DO  111  K=MINK, MAXK, INCRK 

NF1=K-1 

NF2=K*!N-1) 

FINV=FINVER  ! NFl, NF2, ALPHA  ) 

NONCEN=N*STDS 

POW=FPO! ALPHA, NFl, NF2, FINV, NONCEN) 

WRITE! 6,115)  NFl ,NF2, ALPHA, FINV, NONCEN, POW 
115   FORMAT! '0' ,I5,5X,I5,2X,2F10.5,1X,F12. 2,1X,F10.5 ) 

NN!  JJ)=K 

POWER! JJ)=POW 
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JJ=JJ+1 

111   CONTINUE 

PRINT*, 'NN     POWER' 

DO  176  N=1,JJ 

WRITE(6,175)   NN(N),POWER(N) 

175  FORMAT! 'OS  'NN= ' ,F8.3,4X, •POHER= ' ,F10.5) 

176  CONTINUE 

CALL  PLOTTCNN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF((ANS.Eq. 'Y' ).OR.(ANS.EQ. 'Y' ))  GO  TO  500 

GO  TO  999 


C 

C    ANALYZE  ALPHA-LEVEL  VS  POWER  IN  ONE-WAY  ANALYSIS  OF  VARIANCE. 

C 

c 

50     PRINT*, '»  OF  OBSERVATIONS  PER  TREATMENT  (N=  )  ?' 
READ  (5,*)  N 

PRINT*, '«  OF  TREATMENTS  (K=  )?' 
READ  (5,*)  K 

PRINT*,'  STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 
READ  (5,*)  STDS 

PRINT*, 'MAXIMUM   ALPHA   RANGE  ?' 
READ  (5,*)  MAXALP 
PRINT*, 'MINIMUM  ALPHA   RANGE   ?' 
READ  (5,*)  MINALP 

PRINT*, "INCREMENT   ALPHA  VALUE  ?' 
READ  (5,*)  INCRAL 
WRITE(6,295) 
295    F0RMAT(1X,'D0F1   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER') 
JJ=1 
NF1=K-1 
NF2=K*(N-1) 

DO  211  ALPHA=MINALP, MAXALP, INCRAL 
FINV=FINVER  ( NFl ,NF2 , ALPHA ) 
NONCEN=N*ST0S 

P0W=FP0( ALPHA, NFl, NF2,FINV,N0NCEN) 
WRITE( 6,215)  NFl, NF2, ALPHA, FINV,NONCEN,POW 
215   FORMAT! '0' ,I5,3X,I5,2X,2F10 .5,1X,F12. 2,1X,F10.5 ) 
NN( JJ)=ALPHA 
POWER! JJ)=POW 
JJ=JJ+1 
211   CONTINUE 
JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  276  I  =1,JJ 

WRITE! 6,275)   NN! JJ ), POWER! JJ ) 

275  FORMAT!  '0',  'NN= ' ,F8. 3 ,^X, 'POWER= ' ,F10 .5  ) 

276  CONTINUE 

CALL  PLOTTINN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  !Y/N)?' 

READ  !5,7)  ANS 

IF!!ANS.EQ. 'Y' ).OR.IANS.EQ. 'Y' ))  GO  TO  500 


GO  TO  999 


C 

C        ANALYZE  NONCENTRALITY  PARAMETER  VS  POWER  IN  ONE-WAY  ANALYSIS  OF 

C  VARIANCE 

C 

C 

'♦O    PRINT*, '»  OF  OBSERVATIONS  PER  TREATMENT  !N=  )  ?' 

READ  !5,*)  N 

PRINT*, '«  OF  TREATMENTS!  K=)?' 

READ  !5,*)  K 

PRINT*, 'ALPHA-LEVEL   ?' 

READ  15,*)  ALPHA 

PRINTK, 'MAXIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  IS,»)   MAXSD 

PRINT*, 'MINIMUM  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  A 5,*)   MINSD 

PRINTW, 'INCREMENT  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5>»)  INCSD 

WRITE(6,<^95) 
^95    F0RMATI1X,'D0F1   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER') 

JJ=1 

NF1=K-1 

NF2=K*!N-1) 

DO  411  STDS=MINSD,MAXSD, INCSD 

FINV=FINVER  ! NFl ,NF2, ALPHA ) 

NONCEN=N*STDS 

POW=FPO( ALPHA , NFl ,NF2 , FINV ,NONCEN ) 

WRITE! 6,415)  NFl, NF2, ALPHA, FINV, NONCEN,POW 
415   FORMAT! " 0 ' ,IS,3X,I5,2X,2F10.5,1X,F12. 2,1X,F10.S) 

NN! JJ)=NONCEN 
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POWER( JJ)=POW 
JJ=JJ+1 
411  CONTINUE 
JJ=JJ-1 

DO  476  I  =1,JJ 
WRITE(6,475)  NN(I),POWER(I) 

475  FORMAT! '0',  'NN= ' ,F8.3,4X, 'POHER= ' ,F10.5) 

476  CONTINUE 

CALL  PLOTTt NN, POWER, JJ,0) 

PRINT*, "DO  YOU  HANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( (ANS.EQ. 'Y'  ).OR.(ANS.EQ. 'Y'  ))  GO  TO  500 

GO  TO  999 


C 

C        ANALYZE  MULTI-WAY  ANALYSIS  OF  VARIANCE. 

C 

C 


1000  PRINT*, "HOW  MANY  WAYS  ANOVA  CLASSIFICATION?' 

READ(5,*)  M 

DO  600  1=1, M 

PRINT*, 'FOR' , I, '-FACTOR  LEVEL  ,HOW  MANY  LEVELS  ?' 

READ(5,*)  MFII) 
600   CONTINUE 

C 

C 


C        ANALYZE  MODEL  WITH  FACTOR-LEVEL  IN  MULTI-WAY  ANOVA. 

C 

c 

PRINT*, 'IS  THERE  ANY  INTERACTION  ?(Y/N)' 

READ(5,7)  ANS 

IF((ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' ))  GO  TO  3100 

ISUM=0 

IT0T=1 

DO  800  1=1, M 

ISUM=ISUM+(MF(I)-1) 
ITOT=ITOT*MF(I) 
800   CONTINUE 

IERDF=IT0T-ISUM-1 

C 

C 

C        COMPUTE  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT. 

C 

C 

MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  2  RD  FACTOR, PRESS  3)?' 

READ(5,*)  I 

NF1=MF(I)-1 

PICON=FLOAT( ITOT  )/FLOAT( MF( I ) ) 

GO  TO  2000 

C 

C 

C        ANALYZE  MODEL  WITH  FACTOR-LEVEL  AND  2-WAY  INTERACTION  IN 

C  MULTI-WAY  ANOVA. 

C 

C 

5100  PRINT*, 'ARE  THERE  ONLY  2-WAY  INTERACTION  AND  FACTOR-LEVEL?( Y/N ) ' 

READ! 5,7)  ANS 

IF( (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N'  ))  GO  TO  2150 

ISUM1=0 

ISUM2=0 

IT0T=1 

DO  850  1=1, M 

ISUMl=ISUMl+( MF( I)-l ) 
ITOT=ITOT*MF(I) 
850   CONTINUE 

DO  860  1=1, M-1 
DO  870  J=I+1,M 

ISUH2=ISUM2+( MF( I  )-l )«( MFC  J  )-l ) 
870   CONTINUE 
860   CONTINUE 

ISUM=ISUM1+ISUM2 

IERDF=IT0T-ISUM-1 

MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 

PRINT*, 'DO  YOU  WANT  TO  USE  ONLY  FACTOR-LEVEL  (Y/N)?' 

READ(5,7)  ANS 

IF( (ANS.EQ. 'N'  ).0R.( ANS.EQ. 'N' ))  GO  TO  880 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, THEN 
&PRESS  3)?' 

READ(5,*)  I 

NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT( MF( I ) ) 

GO  TO  2000 
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880   PRINT*, 'WHICH  INTERACTION  ANALYSIS?(  IF  YOU  WANT  2  AND  3  WAY  INTER 
FACTION, THEN  PRESS  2   3  ).' 

READ(S,*)  I, J 

NF1  =  (MF(I)-1)*(MF(  J)-l) 

PICON=FLOAT( ITOT )/FLOAT( MF( I  )*MF( J  ) ) 

GO  TO  2000 

C 

C 

C        ANALYZE  MODEL  WITH  FACTOR-LEVEL  AND  3-WAY  INTERACTION  IN 

C  MULTI-WAY  ANOVA. 

C 

C 

2150   PRINT*, 'ARE  THERE  ONLY  3-WAY  INTERACTION  AND  FACTOR-LEVEL  ?(Y/N)' 

READ(S,7)  ANS 

IF( (ANS.EQ. 'N'  ).OR.(ANS.EQ. 'N'  ))  GO  TO  2200 

ISUM1=0 

ISUM3=0 

IT0T=1 

DO  871  1=1, M 

ISUMl=ISUMl  +  ( MF( I  )-l ) 
ITOT=ITOT*MF(I) 

871  CONTINUE 

DO  872  1=1, M-2 

DO  873  J=I+1,M-1 

DO  87"+  K=J+1,M 
ISUM3=ISUM3+( MF( I )-l )*(MF( J )-l )*(  MF(  K )-l ) 

874  CONTINUE 
873   CONTINUE 

872  CONTINUE 
ISUM=ISUM1+ISUM3 
IERDF=IT0T-ISUM-1 
MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 

PRINT*, 'DO  YOU  WANT  TO  ANALYSIS  ONLY  FACTOR-LEVEL  ?(Y/N)' 

READ(5,7)  ANS 

IF(  (ANS.EQ. 'N'  ).0R.( ANS.EQ. 'N'  ))  GO  TO  875 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, PRESS  3)?' 

READ(5,*)  I 

NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT(  MF(  I ) ) 

GO  TO  2000 

875  PRINT*, 'DO  YOU  WANT  TO  ANALYSIS  3-WAY  INTERACTION  ?(Y/N)' 
READ(5,7)  ANS 

IF( (ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' ))  GO  TO  500 

PRINT*, 'WHICH  INTERACTION  ANALYSIS?(  IF  YOU  WANT  2  AND  3  AND  <+  INT 
8ERACTI0N  ,THEN  PRESS  2   3   4  ).' 
READ(5,*)  I,J,K 

NF1  =  (MF(I)-1)*(MF(J)-1)*(MF(K)-1) 
PICON=FLOAT( ITOT )/FLOAT( MF( I  )*MF( J  )*MF( K  ) ) 
GO  TO  2000 


C 

C 

C        ANALYZE  MODEL  WITH  FACTOR-LEVEL  ,2-WAY  INTERACTION  IN  AND 

C  3-WAY  INTERACTION  IN  MULTI-WAY  ANOVA. 

C 

C 

2200  PRINT*, 'ARE  THERE   2-WAY  INTERACTIONS  AND  3-WAY  INTERACTIONS  AND 
SFACTER-LEVEL  (Y/N)?' 
READ(5,7)  ANS 

IF( (ANS.EQ. 'N' ). OR. (ANS.EQ. "N"  ))  GO  TO  999 
ISUM1=0 
ISUM2=0 
ISUM3=0 
IT0T=1 
DO  950  1=1, M 

ISUMl=ISUMl+( MFC  I  )-l ) 
IT0T=IT0T*MF(I) 
950   CONTINUE 

DO  960  I=1,M-1 
DO  970  J=I+1,M 

ISUM2=ISUM2+( MF( I  )-l )»( MFC  J  )-l ) 
970   CONTINUE 
960   CONTINUE 

DO  965  1=1, M-2 

DO  975  J=I+1,M-1 

DO  985  K=J+1,M 
ISUM3=ISUM3  +  ( MFC  I )-l )»C  MFC  J  )-l  )*C  MFC  K )-l ) 
985   CONTINUE 
975   CONTINUE 
965   CONTINUE 

ISUM=ISUM1+ISUM2+ISUM3 

IERDF=IT0T-ISUM-1 

MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT  LEAST  MINIMUM  NUMBER  OF  REPLICATE  PER  TREATMENT  =',MINR 

PRINT*, 'DO  YOU  WANT  TO  ANALYSIS  ONLY  FACTOR-LEVEL  (Y/N)?' 

59 


READ(5,7)  ANS 

IF( (ANS.EQ. 'N'  ).OR.(ANS.EQ. 'N'  ))  GO  TO  990 

PRINT*, 'WHICH  FACTOR  ANALYSIS  (IF  YOU  WANT  3  RD  FACTOR, PRESS  3)?' 

READ(5,»)  I 

NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT( MF( I ) ) 

GO  TO  2000 
990    PRINT*, 'DO  YOU  WANT  TO  ANALYZE  2-WAY  INTERACTION  (Y/N)?' 

READ(5,7)  ANS 

IF( (ANS.EQ.  'N'  ). OR. (ANS.EQ. 'N-  ))  GO  TO  995 

PRINT*, 'WHICH  INTERACTION  ANALYSIS  ?(  IF  YOU  WANT  2  AND  3  WAY  INTE 
&ACTION,THEN  PRESS  2   3  ).' 

READ(S,*)  I, J 

NF1  =  (MF(I  )-!)*(  MF(J)-l) 

PICCN=FLOAT( ITOT )/FLOAT( MF( I )*MF(  J ) ) 

GO  TO  2000 
995    PRINT*, 'DO  YOU  WANT  TO  ANALYZE  3-WAY  INTERACTION  (Y/N)?' 

READ(5,7)  ANS 

IF( (ANS.EQ.  'N'  ).0R.( ANS.EQ.  'N'  ))  GO  TO  500 

PRINT*, 'WHICH  INTERACTION  ANALYSIS?(  IF  YOU  WANT  2  AND  3  AND  ^  INT 
&ERACTION  ,THEN  PRESS  2   3   4  ).' 

READ(5,*)  I,  J,  K 

NF1=(MF(I)-1)*(MF(J)-1)*(MF(K)-1) 

PICON=FLOAT( ITOT )/FLOAT( MF( I )*MF( J )*MF(  K ) ) 

GO  TO  2000 
2000   PRINT*, 'DO  YOU  WANT  TO  PLOT  R(  »  OF  REPLICATE  PER  CELL) 
&VS  POWER  (  Y/N)?' 

READ(5,7)  ANS 

IF( (ANS.EQ.  'Y'  ). OR. (ANS.EQ. 'Y'  ) )  GO  TO  4100 

PRINT*, 'DO  YOU  WANT  TO  PLOT  NONCENTRALITY  VS  POWER  (Y/N)?' 

READ(5,7)  ANS 

IF( (ANS.EQ.  'Y'  ). OR. (ANS.EQ. 'Y' )  )  GO  TO  4200 

PRINT*, 'DO  YOU  WANT  TO  PLOT  ALPHA-LEVEL  VS  POWER  (Y/N)?' 

RE AD (5,7)  ANS 

IF( (ANS.EQ. 'Y'  ). OR. (ANS.EQ. 'Y'  ) )  GO  TO  2300 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( (ANS.EQ. 'Y'  ). OR. (ANS.EQ. 'Y'  ))  GO  TO  500 

GO  TO  999 

C 

C 

C      ANALYZE  NUMBER  OF  SAMPLE  SIZE  PER  TREATMENT  VS  POWER  IN  MULTI-WAY 

C  ANALYSIS  OF  VARIANCE. 

C 

C 

4100   PRINT*, 'ALPHA  LEVEL?' 

READ  (5,*)  ALPHA 

PRINT*, 'STANDARDIZED  SQUARED  DISTANCE  VALUE  ?* 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM  R  VALUE  (  THE  MAXIMUM  VALUE  ON  X-AXIS  )?' 

READ  (5,*)  MAXR 

PRINT*,'  MINIMUM  R  VALUE  (  MUST  BE  MORE  THAN',MINR,'  )?' 

READ  (5,*)  MINR 

PRINT*, 'INCREMENT  R  VALUE  ?' 

READ  (5,*)  INCRR 

WRITE( 6,4095) 
4095   F0RMAT(1X,'D0F1   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER') 

JJ  =  1 

DO  4110  N=MINR, MAXR, INCRR 

NF2=N*IT0T-ISUM-1 

FINV=FINVER  ( NFl ,NF2, ALPHA ) 

NONCEN=N*PICON*STDS 

POW=FPO( ALPHA, NFl ,NF2,FINV,N0NCEN ) 

WRITE( 6,4015)    NFl ,NF2 , ALPHA, FINV,NONCEN,POW 
4015      FORMAK  '0 ' ,I5,3X,I5,2X,2F10.5,1X,F12. 2,1X,F10.5 ) 

NN( JJ)=N 

POWER(JJ)=POW 

JJ=JJ+1 
4110   CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  4076  N=1,JJ 

WRITE (6,4075)   NN(N) ,POWER(N ) 

4075  FORMAT! '0',  'NNs ' ,F8.3,4X, 'POWER= ' ,F10.S ) 

4076  CONTINUE 

CALL  PLOTT(NN,POWER,JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' ))  GO  TO  500 

GO  TO  999 

C 

C 

C     ANALYZE  NONCENTRALITY  PARAMETER  VS  POWER  IN  MULTI-WAY 

C  ANALYSIS  OF  VARIANCE. 

C 
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^200   PRINT*, "HOW  MANY  NUMNER  OF  REPLICATE  (MUST  BE  MORE  THAN',MINR, 
&•  )?■ 

READ(5,»)  IRR 

NF2=IRR»IT0T-ISUM-1 

PRINT*, 'ALPHA-LEVEL   ?' 

READ  (5,*)  ALPHA 

PRINT*, 'MAXIMUM   STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  MAXSD 

PRINT*, "MINIMUM   STANDARDIZED  SQUARED  DISTANCE   VALUE  RANGE  ?' 

READ  (S,*)  MINSD 

PRINT*, 'INCREMENT  STANDARDIZED  SQUARED  DISTANCE  VALUE  RANGE  ?' 

READ  (5,*)  INCSD 

HRITEI 6,4^95) 
4495   FORMATdX, 'DOFl   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER') 

JJ=1 

DO  4911  STDS=MINSD, MAXSD, INCSD 

FINV=FINVER  ( NFl ,NF2 , ALPHA ) 

NONCEN=IRR*PICON*STDS 

POW=FPO(  ALPHA , NFl ,NF2 , FINV ,NONCEN ) 

WRITE! 6,4915)  NFl, NF2 , ALPHA, FINV, NONCEN,POW 
4915   FORMATC  ' 0 ' ,I5,5X,I5,2X,2F10 .5,1X,F12 . 2,1X,F10.S  ) 

NN( JJ)=NONCEN 

POWER( JJ)=POW 

JJ=JJ+1 
4911   CONTINUE 

JJ=JJ-1 

DO  1476  I  =1,JJ 

WRITE( 6,1475)   NN( I ) ,POWER( I ) 

1475  FORMAT!  '0',  'NN= ' ,F8. 3 ,4X, 'POWER= ' ,F10 .5 ) 

1476  CONTINUE 

CALL  PLOTTINN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,7)  ANS 

IF((ANS.EQ. 'Y' ).OR.(ANS.EQ. 'Y' ))  GO  TO  500 

GO  TO  999 

C 

C 

C     ANALYZE  ALPHA-LEVEL  PARAMETER  VS  POWER  IN  MULTI-WAY 

C  ANALYSIS  OF  VARIANCE. 

C 

c 

2300   PRINT*, 'HOW  MANY  NUMNER  OF  REPLICATE  (MUST  BE  MORE  THAN',MINR, 
&'   )?' 

READ(5,*)  IRR 

NF2=IRR*IT0T-ISUM-1 

READ  (5,*)  K 

PRINT*, 'STANDARDIZED  SQUARED  DISTANCE  VALUE  ?' 

READ  (5,*)  STDS 

PRINT*, 'MAXIMUM   ALPHA   RANGE  ?' 

READ  (5,*)  MAXALP 

PRINT*, 'MINIMUM  ALPHA   RANGE   ?' 

READ  (5,*)  MINALP 

PRINT*, 'INCREMENT   ALPHA  VALUE  ?' 

READ  (5,*)  INCRAL 

WRITE( 6,1295) 
1295   FORMATdX, 'DOFl   D0F2   ALPHA-LEVEL   F-INVERSE   NONCENTRAL   POWER*) 

JJ=1 

NF1=K-1 

NF2=K*(N-1) 

DO  1211  ALPHA=MINALP, MAXALP, INCRAL 

FINV=FINVER  ( NFl ,NF2, ALPHA  ) 

NONCEN=IRR*PICON*STDS 

POW=FPO( ALPHA, NFl, NF2, FINV, NONCEN) 

WRITE( 6,1215)  NFl, NF2, ALPHA, FINV, NONCEN, POW 
1215   FORMAT! '0 ' ,IS,3X,I5,2X,2F10.5,1X,F12. 2,1X,F10.5 ) 

NN(JJ)=ALPHA 

POWER(JJ)=POW 

JJ=JJ*1 
1211      COTfTINUE 

JJ=JJ-1 

PRINT*) "NNCJJ )■ ,MINK,MAXK,INCRK 

DO  1276  I  =1,JJ 

WRITE(6,1275)   NN( JJ  ), POWER! JJ  ) 

1275  FORMAT! 'C,  ■NN=*  ,F8.3,4X, 'POWERs' ,F10.5) 

1276  CONTINUE 

CALL  PLOTT(NN,POWER,JJ,0) 
PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 
READ  (5,7)  ANS 

IF((ANS.EQ. 'Y' J.GR.CANS.EQ. 'Y* ))  (K)  TO  500 
GO  TO  999 
999   STOP 
END 

* 
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»  APPROXIMATION   TO  THE  F-  INVERSE  DISTRIBUTION.  » 

*  CALCULATE  F-INVERSE  GIVEN  DOFl ,D0F2, ALPHA.  » 

*  Q(  F-ALPHA  |D0F1,D0F2  )=  ALPHA.  * 

*  * 

FUNCTION  FINVER(N1,N2,ALP) 
REAL  ALP,NUMER,DENUM,LAM,XP,YP,A,B,FIN1 
DATA  CO ,C1 ,C2/2 . 515517 ,0 . 802852 ,0 . 010328/ 
DATA  Dl, 02, D5/1. '♦32788, 0.189269, 0.001308/ 
B=Nl/2. 
A=N2/2. 
C 

c 

C     CALCULATE  NORMAL  QUANTILE. 

C 

C 

T=SQRT(AL0G(1./(ALP**2.  ))) 

NUMER=C0+C1*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2+D3*T»»3 

XP=T-(NUMER/DENUM) 


C 

C     CALCULATE  THE  INVERSE  F-DISTRIBUTION. 

C 

C 


LAM=(XP**2-3. )/6. 

W2=  (l./(2.*B-l.  ))-(l./(2.*A-l.  )) 

H^=    <  l./(2.*A-l.  ))  +  (l./(2.*B-l.  )) 

H  =  2.*W'***(-1.  ) 

Wl={LAM+(5./6.  )-(2./(3.*H))) 

W3=(XP*(H+LAM)**.S)/H 

W=(W3-W2*H1)*2 

IF(W.LT.174)  THEN 
FIN1=EXP(H) 

ELSE 

FIN1=1000 

ENDIF 
FINVER=FIN1 
RETURN 
END 

^  ^  KK  K.^  ^  A  A  W  ^  A  A  A  ^AA^^AWAW^WAA^  AAAWAAAAAAAATWKRWKWWAXWWWKKWKWWA^ 

*  * 

*  GIVEN  DOFl, D0F2,FINV, NON-CENTRAL  PARAMETER,  * 

*  THEN  COMPUTE  POWER.  * 

*  (I.E.   CALCULATE  NON-CENTRAL  F-DISTRIBUTION  (  C.D.F.)   * 

*  * 

AAAAAAWWK^WJWAWKWW^WWWAIWWW^^AA  AAAA  AAAA  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

FUNCTION    FP0(ALP,N1,N2,FIN,N0N) 
REAL  NON 
INTEGER  N1,N2 

DATA  A1,A2, A3, P/. ^361836, -.1201676,. 9372980,. 33267/ 
CONST   =  1.0  /  SQRT  (2.0  *  3.1^15927) 


C 

C 

C     CALCULATE  THE  NORMAL  APPROXIMATION 

C 

C 


C1=(N1*FIN)/(N1+N0N) 

C2=l-(2./(9.*N2)) 

C3=(2.*(N1+2.*N0N))/(9.*(N1+N0N)*«2) 

CCl  =  (Cl**(l./3.  ))*C2-(1-C3) 

CC2  =  (  (C3  +  (2./(9.*N2))*(Cl)**(2./3.  ))»».5) 

FX=CC1/CC2 

C 

C 

C     CALCULATE  THE  NONCENTRAL  F-DISTRIBUTION  USING  NORMAL  CDF. 

C 

C 

COF  B  0.0 

IF(FX  .EQ.O)  THEN 

CDF  =  0.5 

ELSE  IF(FX.GT.O)  THEN 

T=1./(1.+P*FX) 

XSQU=(FX»FX)*.5 

IF(XSQU.LT.17't)  THEN 

ZPDF=CONST»EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 
C  ZPDF=CONST*EXP( -XSQU ) 

CDF=1-ZPDF*( A1»T+A2»T»»2+A3»T»»3 ) 

ELSE 

FX=-FX 

T=1./(1.+P*FX) 
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XSQU=(FX»FX)*.5 

IF(XSQU.LT.17<i)  THEN 

ZPDF=CONST*EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF»( A1*T+A2*T»*2+A3*T**3 ) 
END  IF 
FP0=1-CDF 
RETURN 
END 

*  * 

*  GIVEN  ITOT  AND  ISUM  ,       THEN   CALCULATE  MINIMUM  REPLICATE.     * 

FUNCTION  MINIR(ITOT,ISUM) 
INTEGER  IR,  ITOT , ISUM, ICHECK ,MINIR 
IR=2 
10    ICHECK=IR*IT0T-ISUM-1 
IF( ICHECK. LE.O)  THEN 
IR=IR+1 
GO  TO  10 
ELSE 

MINIR=IR 
ENDIF 
RETURN 
END 
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APPENDIX  G 
THE  POWER  OF  T-TEST  PROGRAM  LIST 


WW  WWTTTTTt  ^  K  R  K^  K  KK  ft  KK^^  WW  WWWWWWWW  W^rWW  W  WITW  W  WW  W  W  WW  WW^W  AAAAAAAAXAKXWWWWWWWAJIJV 


« 

« 

* 
* 
* 

POWER  ALGORITHM  (  T  TEST  ) 

* 

* 

DIRECTED  BY  :  PROFESSOR   DONALD. R. 

BARR 

WRITTEN  BY  :  HUR,  SEONG  PIL 

JULY 

1986 

* 

DEPARTMENT  OF   O.R. 

NAVAL 

POSTGRADUATE 

SCHOOL 

* 

PROGRAM  IS  USED  THE  PRESENTATION 

OF  POWER 

ALGORITHM  OF 

T-TEST. 

* 

* 

(1)  ONE-SAMPLE  CASE  TEST. 

* 

* 

.  ONE-SIDED  TEST  (  HA  : 

MU>MUO  ) 

* 

a 

.  ONE-SIDED  TEST  (  HA  : 

MU<MUO  ) 

w 

a 

.  TWO-SIDED  TEST  (  HA  : 

MU  MUO  ) 

* 

« 

(2)  TWO-SAMPLE  CASE  TEST. 

« 

* 

.  ONE-SIDED  TEST  (  HA  : 

MU>MUO  ) 

* 

* 

.  ONE-SIDED  TEST  (  HA  : 

MU<MUO  ) 

* 

* 

.  TWO-SIDED  TEST  (  HA  : 

MU  MUO  ) 

* 

* 

NOTE  :  Subroutine  PLOTT  is  NONIMSL 

subroutine  library. 

w 

» 

* 

AAA^RA^AARAAA  A  AWARAAAA^AWRWW  WWAWXWWWA  RWRWWXWWWWWWWAXARRWKWWWWAAARA  JItWWWT 

REAL  NN( 200  )  ,POWER( 200  )  ,DF ,NONCDF ,TCDF1 ,TCDF2 ,TINV1 
REAL  TINV,NONCEN,TX, ALPHA, CDF, CONST  ,MAXMU1 ,MINMU1 ,INCMU1 
REAL  MAXDEL,MINDEL,INCDEL,P0W,P0W1,P0W2,MUZER0,STDV,MU1 
INTEGER  JJ,K,N 
CHARACTER*1  ANS 

C 

C 

C 

C 

C 

100 

10 


INPUT  USER  OPTIONS. 


PRINT*, 'DO  YOU  WANT  TO  TEST  ONE  SAMPLE  T-TEST( Y/N)? ' 

READ(5,10)  ANS 

FORMAT  (AD 

IF( (ANS.EQ. 'N' ).OR.(ANS.EQ. 'N' ))  GO  TO  160 

PRINT*,'  SAMPLE  SIZE  ?• 

READ(5,*)  N 

PRINT*,'  ALPHA-LEVEL  ?' 

READ(5,*)  ALPHA 

PRINT*,'  MU-ZERO  VALUE   ?' 

READ(5,*)  MUZERO 

PRINT*,'  STANDARD  DEVIATION  VALUE  ?' 

READ(5,*)  STDV 


ONE -SAMPLE  ONE-SIDED  TEST,  HA  :  MU  >  MUO. 


'N'  ))  GO  TO  30 
HO:MU=MUO  VS  HA: 


MU>MUO  (Y/N)?' 


PRINT*, 'DO  YOU  WANT  TO  USE  ONE-SIDED  TEST  (Y/N)?' 
READ! 5,10)  ANS 
IF( (ANS.EQ. 'N' ). OR.  (ANS.EQ, 
PRINT*, 'DO  YOU  WANT  TO  TEST 
READ(5,10)  ANS 
IF( (ANS.EQ. 'N* ). OR. (ANS.EQ. 
PRINT*,'  MAXIMUM  t*Jl  VALUE  i 
READ  (5>»)  MAXMUl 

PRIhTT*,'  MINIMUM  MUl  VALUE  (  THE  MINIMUM  VALUE  OF  X-AXIS  )? 
&  (COra)ITION  :MU1  MUST  BE  MORE  THAN  MU-ZERO  )' 
READ  (5>«)  MINMUl 

PRINT*,'  LARGE  INCREMENT   N  VALUE  ?' 
READ  (5>M)  INCMUl 
N=N-1 


N'  )) 
THE 


GO  TO  '+0 
MAXIMUM  VALUE 


OF  X-AXIS  )?' 


GENERATE  T-INVERSE  VALUE  (  CRITICAL  REGION  VALUE  ). 


TINV=TINVER(N, ALPHA) 
JJ=1 
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DO  50  MU1=MINMU1,MAXMU1,INCMU1 

C 

C 

C     COMPUTE  NON-CENTRALITY  VALUE. 

C 

C 


DF=N 

NONCEN=( ( MUl-MUZERO  )*SQRT( DF  +  1 )  )/STDV 


C 

C      GENERATE  CDF  OF  NONCENTRAL  T  DISTRIBUTION. 

C 

C 


TCDF=NONCDF( N,NONCEN,TINV ) 

P0W=1-TCDF 

NN( JJ)=MU1 

POWER! JJ)=POW 

JJ=JJ+1 
50     CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  60  1=1, JJ 

WRITE(6,70)   NN(I),POWER(I) 
70    FORMAT! '0',  'NN= ' ,F8.3 ,4X, ' POWER= ' ,F10.5 ) 
60    CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF((ANS.Eq. 'Y'  ).OR.(ANS.EQ. 'Y' ))  GO  TO  100 

GO  TO  80 


C 

C     ONE -SAMPLE  ONE-SIDED  TEST,  HA  :  MU  <  MUO. 

C 

C 


'♦O    PRINT*, 'DO  YOU  WANT  TO  TEST  HO:MU=MUO  VS  HA:MU<MU0  (Y/N)?' 

READ(5,10)  ANS 

IF( (ANS.EQ. 'N'  ).OR.(ANS.EQ. 'N'  ))  GO  TO  80 

PRINT*,'  MAXIMUM  MUl  VALUE  (  THE  MAXIMUM  VALUE  OF  X-AXIS  )?' 

READ  (5,*)  MAXMUl 

PRINT*,'  MINIMUM  MUl  VALUE  (  THE  MINIMUM  VALUE  OF  X-AXIS  )? 
&  (CONDITION  :MU1  MUST  BE  LESS  THAN  MU-ZERO  )  ' 

READ  (5,*)  MINMUl 

PRINT*,'  INCREMENT   MUl  VALUE  ?' 

READ  (5,*)INCMU1 

N=N-1 

TINV=TINVER(N, ALPHA) 

TINV=-TINV 

JJ  =  1 

DO  90  MU1=MINMU1, MAXMUl, INCMUl 

DF=N 

NONCEN=( ( MUl-MUZERO  )*SQRT( DF+1 )  )/STDV 

TCDF=NONC0F( N,NONCEN,TINV ) 

POW=TCDF 

NN( JJ)=MU1 

POWER! JJ)=POW 

JJ=JJ+1 
90     CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  110  I=1,JJ 

WRITE(6,120)   NN(I),POWER(I) 
120    FORMAT!  '0',  'NN= ' ,F8.3,4X, 'POWER= ' ,F10.S ) 
110   CONTINUE 

CALL  PLOTT!NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF! (ANS.EQ. 'Y' ) .OR. (ANS.EQ. 'Y'  ))  GO  TO  100 

GO  TO  80 

C 

C 

C     ONE-SAMPLE  TWO-SIDED  TEST,  HA  :  MU  =/  MUO. 


30   PRINT*, 'DO  YOU  WANT  TO  USE  TWO-SIDED  TEST! Y/N)?' 
READ. 5,10)  ANS 

IF! (ANS.EQ. 'N'  ). OR. ! ANS.EQ. 'N'  ))  GO  TO  80 

PRINT*, 'MAXIMUM  MUl  VALUE  !  THE  MAXIMUM  VALUE  ON  X-AXIS  )?' 
READ  (5,*)  MAXMUl 

PRINT*,'  MINIMUM  MUl  VALUE  !  THE  MINIMUM  VALUE  ON  X-AXIS  )? 
S  (CONDITION  :MU1  MUST  BE  LESS  THAN  MU-ZERO  )' 
READ  (5,*)  MINMUl 
PRINT*,'  INCREMENT   MUl  VALUE  ?' 
READ  (5,*)INCMU1 
N=N-1 
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ALPHA=ALPHA/2. 

TINV=TINVER( N, ALPHA ) 

JJ  =  1 

DO  120  MU1=MINMU1,MAXMU1,INCMU1 

DF=N 

NONCEN=((MUl-MUZERO)itSQRT(DF  +  l.  )  )/STDV 

TCDF1=N0NCDF(N,N0NCEN,TINV) 

P0H1=1.-TCDF1 

TINV1=-TINV 

TCDF2=N0NCDF(N,N0NCEN,TIN'v'l) 

POW2=TCDF2 

P0W=P0H1+P0H2 

NN( JJ)=MU1 

POWER( JJ)=POW 

JJ=JJ+1 
130    CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  1^0  1  =  1, JJ 

WRITE( 6,150)   NN(I),POWER(I) 
150    FORMAT! '0',  ■NN= ' ,F8. 3 ,4X, ' POWER= ' ,F10 .5 ) 
1^0   CONTINUE 

CALL  PLOTT(NN,POHER,JJ,0) 

PRINT*, '00  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF  I (ANS.EQ. 'Y'  ).OR.(ANS.EQ. 'Y'  ))  GO  TO  100 

GO  TO  80 

C 

C 

C     TWO-SAMPLE  ONE-SIDED  TEST,  HA  :  MUX  >  MUY. 

C 

c 


160   PRINT*, 'DO  YOU  WANT  TO  TEST  TWO  SAMPLE  T-TEST( Y/N)? ' 

READ(5,10)  ANS 

IF( (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N'  ))  GO  TO  80 

PRINT*,'  SAMPLE  SIZE  FROM  POPULATION  1  ?' 

READ(5,*)  N 

PRINT*,'  SAMPLE  SIZE  FROM  POPULATION  2  ?' 

READ(5,*)  M 

PRINT*,'  ALPHA-LEVEL  ?' 

READ(5,*)  ALPHA 

PRINT*, 'DO  YOU  WANT  TO  USE  ONE-SIDED  TEST  (Y/N)?' 

READ(5,10)  ANS 

IF( (ANS.EQ.  'N'  ). OR. (ANS.EQ. 'N'  ))  GO  TO  170 

PRINT*, 'DO  YOU  WANT  TO  TEST  HO:  MUX=MUY  VS  HA:  MUX>MUY  (Y/N)?' 

RE AD (5,10)  ANS 

IF( (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N' ))  GO  TO  180 

PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 

READ  (5,*)  MAXDEL 

PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 
C     8  (CONDITION  :MU1  MUST  BE  MORE  THAN  MU-ZERO  )?' 

READ  (5,*)  MINOEL 

PRINT*,'  INCREMENT  NONCENTRALITY  VALUE?' 

READ  (5,*)INCDEL 

Nl=M+N-2 

TINV=TINVER( Nl , ALPHA ) 

JJ  =  1 

DO  190  NONCEN=MINDEL, MAXDEL, INCDEL 

TCDF=NONCDF( Nl ,NONCEN,TINV ) 

P0W=1-TCDF 

NN( JJ)=NONCEN 

POWER( JJ)=POW 

JJ=JJ+1 
190   CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  200  1  =  1, JJ 

WRITE( 6,210)  NN(I),POWER(I) 
210   FORMAK'C,  '1^=' ,F8.3,4X, 'POWER=' ,F10.5) 
200   CONTINUE 

CALL   PLOTKNN, POWER, JJ,0) 

PRINT»,*DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5>10)  ANS 

IF( (ANS.EQ. 'Y'  ). OR. (ANS.EQ. 'Y'  ))  GO  TO  100 

GO  TO  80 

C 

C 


C     TWO-SAMPLE  ONE-SIDED  TEST,  HA  :  MUX  <  MUY. 
C 

c 


180   PRINT*, 'DO  YOU  WANT  TO  TEST  HO:MUX=MUY  VS  HA:MUX<MUY  (Y/N)?' 
RE AD (5,10)  ANS 

IF( (ANS.EQ.  'N'  ). OR. (ANS.EQ. 'N'  )  )  GO  TO  80 
PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 
READ  (5,*)  MAXDEL 
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PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 

READ  (5,*)  MINDEL 

PRINT*,'  INCREMENT  NONCENTRALITY  VALUE?' 

READ  (5,*)INCDEL 

Nl=M+N-2 

TINV=TINVER( Nl , ALPHA ) 

TINV2=-TINV 

JJ  =  1 

DO  220  NONCEN=MINDEL,MAXDEL,INCDEL 

TCDF=NONCDF( Nl ,N0NCEN,TINV2 ) 

P0W=TCDF 

NN( JJ)=N0NCEN 

PGWER( JJ)=POH 

JJ=JJ+1 
220   CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  230  1=1, JJ 

WRITE(6,240)   NN(I  ),POWER(I) 
240    FORMAT('0',  'NN= ' ,F8. 5 ,4X, ■ POWER= ' ,F10.5 ) 
250   CONTINUE 

CALL  PLOTT(NN,POHER,JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF( (ANS.EQ.  'Y'  ).OR.(ANS.EQ. 'Y'  ))  GO  TO  100 

GO  TO  80 


C 

C     TWO-SAMPLE  TWO-SIDED  TEST,  HA  :  MUX  =/  MUY. 
C 

c 


170   PRINT*, 'DO  YOU  WANT  TO  USE  TWO-SIDED  TEST (Y/N)?' 

READ(5,10)  ANS 

IF( (ANS.EQ. 'N'  ). OR. (ANS.EQ. 'N' ))  GO  TO  80 

PRINT*,'  MAXIMUM  NONCENTRALITY  VALUE  ?  ' 

READ  (5,*)  MAXDEL 

PRINT*,'  MINIMUM  NONCENTRALITY  VALUE?' 

READ  (5,*)  MINDEL 

PRINT*, 'INCREMENT   NONCENTRALITY  VALUE?' 

READ  (5,*)INCDEL 

Nl=M+N-2 

ALPHAl=ALPHA/2. 

TINV=TINVER( Nl , ALPHA! ) 

JJ=1 

DO  250  NONCEN=MINDEL, MAXDEL, INCDEL 

TCDF1=N0NCDF(N1,N0NCEN,TINV) 

P0W1=1.-TCDF1 

TINV1=-TINV 

TCDF2=N0NCDF(N1,N0NCEN,TINV1) 

POW2=TCDF2 

P0W=P0W1+P0W2 

NN( JJ)=NONCEN 

POWER( JJ)=POW 

JJ=JJ+1 
250   CONTINUE 

JJ=JJ-1 

PRINT*, 'NN       POWER' 

DO  260  1=1, JJ 

WRITE(6,270)   NN( I ) ,POWER( I) 
270   FORMAT('0',  'NN=' ,F8.3,4X, 'POWER=' ,F10.5) 
260   CONTINUE 

CALL  PLOTT(NN, POWER, JJ,0) 

PRINT*, 'DO  YOU  WANT  TO  PLOT  AGAIN  (Y/N)?' 

READ  (5,10)  ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' ))  GO  TO  100 

GO  TO  80 
80   STOP 

END 

*  APPROXIMATION  TO  THE  T-  INVERSE  DISTRIBUTION. 

»  CALCULATE  T-INVERSE  GIVEN  DOF  ALPHA. 

»  Q(  T-ALPHA  I  DOF  )=  ALPHA. 
* 

FUNCTION  TINVER(N, ALPHA) 

REAL  ALPHA ,NUMER ,DENUM,LAM ,XP ,G1 ,G2 ,G3 ,X,TINVER 

INTEGER  N 

DATA  CO ,C1 ,C2/2 .515517,0 .802853 ,0 . 010328/ 

DATA  Dl ,D2 ,D3/1 . 432788 ,0 . 189269 ,0 . 001308/ 


C_ 

c 

C     CALCULATE  NORMAL  QUANTILE. 

C 

C 
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T=SQRT(AL0G(1./(ALPHA*»2.  ))) 

NUMEk=C0+C1*T+C2*T**2 

DENUM=1 . +01»T+D2*T**2 . +D3*T**3 . 

XP=T-(NUMER/DENUM) 

C 

C 

C     CALCULATE  THE  INVERSE  T-DISTRIBUTION. 

C 

C 


X=XP 

G1=.25»(X»*3+X) 

G2  =  ( 1./96.  )*(5.*X**5+16.*X**5+3.*X) 
G3  =  ( l./38^.  )*(3.*X*i*7+19.*X**3  +  17.*X**3-15.*X) 

G'+=(  1 .  /92160  .  )*(  79 .  »(X**9+776  .  *X**7+lCh82  .  *X**5-1592  .  *X**3-9<+5 .  *X  ) 
TP=X  +  (  Gl/N  )  +  (  G2/(  N*N  ) )  +  (  G3/(  N**3  ) )  +  (  64/ (  N**"^  ) ) 
TINVER=TP 
RETURN 
END 
iHi*iHt**-)i*****mHi*******      INFORMATION  ***it*************x*)()(*xxnxxx')in**** 

*  THE  OBJECTIVE  OF  THIS  PROGRAM  IS  TO  CALCULATE  NON-CENTRAL  T    * 
»         DISTRIBUTION  (CDF)  * 

*  THE  CALCULATION  METHOD  IS  NORMAL  APPROXIMATION.  * 

*******»*»*»*i*********  VARIABLE  DECLARATION  ***jf»»*****x)0(it»>t»»»xxxKX* 

FUNCTION  NONCDF(NF,TLAM,TSTAR) 
INTEGER  NF 

REAL  TSTAR,TLAM,TX,  NONCDF,CDF 

DATA  Al , A2 ,A3 ,P/ . '+361836 , - . 1201676 , . 9372980 , , 33267/ 
CONST   =  1.0  /  SQRT  (2.0  *  3.1415927) 


C 

C     CALCULATE  THE  NORMAL  APPROXIMATION 

C 

C 


Cl=( l.-( l./(4.*NF  ) )  )*TSTAR-TLAM 

C2=SQRT( 1 . +( TSTAR**2/( 2 .*NF  )  ) ) 

TX=C1/C2 

CDF  =  0.0 

IF(TX  .EQ.O)  THEN 

CDF  =  0.5 

ELSE  IF(TX.GT.O)  THEN 

T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 

IF(XSQU.LT.17'+)  THEN 

ZPDF=CONST*EXP( -XSQU) 

ELSE 

ZPDF=0 

ENDIF 

CDF=1-ZPDF*( A1*T+A2*T»*2+A3»T**3 ) 

ELSE 

TX=-TX 

T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 

IF( XSQU. LT. 174)  THEN 

2PDF=C0NST»EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF*( A1*T+A2»T**2+A3*T**3 ) 
END  IF 

NONCDF=CDF 
RETURN 
END 
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